Mathematics of Computation 


A Quarterly Journal 
Edited by 
ALAN FLETCHER P. C. HAMMER EUGENE ISAACSON 
Y. L. LUKE DANIEL SHANKS A. H. TAUB 
R. 8S. VARGA J. W. WRENCH, JR. D. M. YOUNG 


HARRY POLACHEK, Chairman 


XVI 
NOS. 77-80 
1962 


Formerly: Mathematical Tables and other Aids to Computation 


Published for the 
National Academy of Sciences—National Research Council 
By the 


American Mathematical Society 











Mathematics 





q/: 
47 





mC 
the a; 
N 


| 
| 
| 
| 





it a et aa 


~ 7 ina bh 
ae" 


aes S ; $ 
+7 
+ U ww HIGAN 


JAN 1i 1962 


Mathematics “= 
of 
Computation 





A journal devoted to advances in numerical analysis, 
the application of computational methods, mathematical tables, 
high-speed calculators and other aids to computation 





Formerly: Mathematical Tables and other Aids to Computation 


Published Quarterly for the 
National Academy of Sciences—National Research Council 
By the 
American Mathematical Society 





Editorial Committee 
Division of Mathematics 
National Academy of Sciences—National Research Council 
Washington, D. C. 


H. Fd sag nev ome, _ Mathematics Laboratory, David Taylor Model 
Basin, Washington 7, D 
ALAN FLETCHER, Univcrsity of Liverpool, Liverpool 3, England 
P. C Hammer, University of Wisconsin, Madison 6, Wisconsin 
Evucens Isaacson, New York University, New York 3, New York 
Y. L. Luxg, Midwest Research Institute, Kansas City ‘10, Missouri 
Dante. SHanxs, Applied Mathematics Laboratory, David Taylor Model Basin, 
Washington 7, D. C. 
A. H. Tavs, University of Illinois, Urbana, Illinois 
R. 8. Varea, Case Institute of Technology, Cleveland 6, Ohio 
J. W. Wrencu, Jr., Applied Mathematics Laboratory, David Taylor Model Basin, 
Washington 7, D.C. 
D. M. Youne, Computing Center, University of Texas, Austin 12, Texas 


Information to Subscribers 


The journal is published quarterly in one volume per year with issues numbered 
seriaily since Volume I, Number 1. Starting with January, 1959 subscriptions are 
$8.00 per year, single copies $2.50. Other issues are available as follows: 
Volume I (1943-1945), Nos. 10 and 12 only are available; $1.00 per issue. 
Volume II (1946-1947), Nos. 13, 14, 17, 18, 19, and 20 only available; $1.00 


per issue. 
Volume IIT (1948-1949), Nos. 21-28 available. $4.00 per year (four issues), 
$1.25 per issue. 
Volume IV—XII (1950 through 1958), all issues available; $5.00 per year, $1.50 
per issue. 
Microcard Edition 
Volumes I—X (1943-1956), Nos. 1-56 are now available on Microcards and may be 


purchased from the Microcard Foundation, Box 2145, Madison 5, Wisconsin, at a 
cost of $20.00 for the complete set. Succeeding volumes are available on request. 


Information to Contributors 


All contributions intended for publication in Mathematics of Computation and all 
books for review should be addressed to H. Polachek, Technical Director, Applied 
Mathematics Laboratory, David Taylor Model Basin, Washington 7, D. C. The 
author may suggest an appropriate editor for his paper. Manuscripts should be 
typewritten double-spaced in the format used by the journal. For journal abbrevia- 
tions, see Mathematical Reviews, v.21, Index, 1960. Authors should submit the orig- 
inal and one copy, and should retain one copy. 


Subscriptions, address changes, business communications and payments should be 
sent to: 


AMERICAN MATHEMATICAL SOCIETY 

190 Hope Street 

Providence 6, Rhode Island 

Published Quarterly for the 
NATIONAL ACADEMY OF SCIENCES—NATIONAL RESEARCH COUNCIL 
By the 
AMERICAN MATHEMATICAL SOCIETY 
Baltimore, Maryland and Providence, Rhode Island 


Copyright © 1962 by the American Mathematical Society 
Printed in the United States of America 
Second-class postage paid at Baltimore, Md. 





o 





SY 
2260 
BF Y/K 
Correlations and Spectra for Non-Stationary 
Random Functions 


By J. Kampé de Fériet and Francois N. Frenkiel 


1. Introduction. The definition of the correlation function and the spectrum of a 
stationary random function is now classical, but in many applications one feels the 
need to extend this definition to random functions, which, although non-stationary, 
are in some sense nearly stationary. We suggest, therefore, for the definition of the 
correlation of a random function whose covariance I(t, s) is known, the limit 


{Al 


a oe we. i h 
(1.1) R(h) = lim 7 he r(é _ > & + *) dé, 


T++0 


if this limit exists for every h. The spectrum S(A) can then be obtained from R(h) 
in the classical way. 

We are led to the above definition of the correlation function R(h) by the follow- 
ing considerations: we determine the sample-correlation from a truncated sample of 
the random function; we then obtain a sub-correlation, Rr(h), of the random func- 
tion (defined as the correlation of the truncated random function) by averaging 
the sample correlations; finally, the correlation R(h) is defined by (1.1) as the 
limit of Rr(h), if this limit exists. 

The function R(h), so defined, has all the properties of a correlation function. 
If the random function is stationary (wide sense) [4, p. 95-96], our definition coin- 
cides with the classical definition. The estimation of the correlation of a stationary 
random function has been considered extensively in the literature, particularly by 
U. Grenander and M. Rosenblatt [6], R. B. Blackman and J. W. Tukey [1], and E. 
Parzen [13, 14]. In order to evaluate how good the estimate R(h) is from the sample- 
correlations pr(h, w), which are the only experimental observables, we compute the 
variance of the random variables pr(h, w) about Rr(h), and then we compute (fora 
fixed h) an upper bound of R(h) — Rr(h) for large T. 

This paper is especially concerned with the case in which the random function 
has a periodic covariance ['(t + 7, s + 7) = I(t, s). To appreciate the scope of 
the above condition, let us note that it is always satisfied when the random function 
is a sum of two uncorrelated random functions, one being a stationary (wide sense ) 
random function and the other a periodic random function. 

The last part of the paper is devoted to the estimate of R(h) for a non-stationary 
random step-function V(t, w), similar to the one introduced by N. Wiener, 


(1.2) Vit, w) = X.(w), n—-1St<n,n=1,2,---, 


where [Xi(w), --- X,(w), ---] is a sequence of independent random variables 
taking only the values —1 and +1 with equal probability. 

An experimental function will be constructed using a table of random numbers, 
and sample-correlations will be determined. Estimates for the sub-correlations are 
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then determined by taking averages over the experimental correlations. The accu- 
racy of these estimates is characterized by giving the variance for the departure of 
the sample-correlations from the estimated sub-correlations. The experimental data 
are then compared with the theoretical results. (Cf. [5], [9].) 


2. Random Functions, Covariances and Correlation Functions. We consider here 
a random function of a real variable t as an ensemble of real functions f(t, w), where 
w is a parameter chosen at random in some set 2 according to a probability measure 
u [7]. A sample of the random function f(t, w) is simply the real function f(t, wo) 
corresponding to a particular choice of w in the set Q. It is convenient for many 
applications to take for 2 a function space; each point w is then a function w(t) 
belonging to some prescribed class of functions (e.g., a continuous function on 
{0, 1]). One has thus for each sample f(t, #) = w(t). When this particular choice is 
made for Q, one says that the random function is of ‘function space type’’ [4, p. 67]. 

The following general hypotheses shall be made with regard to the random 
functions considered in the present paper: 

H, . f(t, w) is measurable with respect to the product measure m X yu (where m 
is the Lebesgue measure on the real line —~ <t < +). 

H. . For each t, f(t, w) € L’(Q): f(t, w)? < +. (If F(w) € L(Q) we denote its 


mean value by [ F(w)du = F(w).) 
Q 





H; . H: implies f(t, w) € L(Q); we suppose f(t, w) = 0. 
H, . It follows from Hy, that the covariance 


(2.1) I(t, s) = f(t, w)f(s, w) 
exists for all ¢’s and s’s. We assume that, 
(2.2) T(t, t) € Lfa, 6] 


for every finite interval a < ¢ < b. From (2.2), by the Fubini-Tonelli theorem, 
[8, Vol. 1, p. 609] it follows that 


(2.3) f(t, #) € L*{a, bj 


for almost all samples, in any finite interval a < t < b. This implies also that 


(2.4) f(t, w) € Lfa, 6] for almost all samples. 
In our earlier paper [9] instead of (2.2), we assumed that 
(2.2’) F(t, s) € L(d’) 


for every finite rectangle A’ in the plane (t, s). We are indebted to the referee for a 
simple counter example showing that (2.2’) does not always imply (2.3) and for 
the suggestion that we replace (2.2’) by (2.2); the proof that (2.2) implies (2.3) by 
the Fubini-Tonelli theorem is straightforward. 

If the random function f(t, w) is stationary 


(2.5) I(t, s) = p(t — s) 


where p(h) is called a correlation function. Due to H, , p(h) is uniformly continuous 
in any finite interval [3]. A real function p(A) is a correlation function if, and only 
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if, it is symmetric and positive-definite 


(2.6) p(—h) = ph) 

(2.7) DD X;Xe(hj — he) 2 0 
3 

for any set [h; , --- h,] with n arbitrary. 


3. Truncated Samples of Random Functions. In experiments concerning a ran- 
dom function f(t, #) one can materialize, as a rule, the sample of the function only for 
a finite interval, that is, one knows only truncated samples. As far as finite intervals 
are considered, one often uses the notation [— 7’, +7'] for the interval in which the 
samples are known in the experiment. Rather than this two-sided (symmetric with 
respect to t = 0) truncation, we shall prefer here a one-sided truncation (starting at 
t = 0) and we will define a truncated sample by 
Sr(t, wo) = f(t, wo), Osis Ze 
(3.1) 

fr(t, wo) = 0, t<Oort> T. 


This definition implies that the experiment starts at t = 0; we assume that it could 
be extended for an arbitrary time T in the future, but not in the past (time prior 
to the beginning of the experiment). From the samples we will draw some inference 
with regard to the random function for 0 < t < +, but completely ignore it for 
t< 0. 


4. Correlation and Spectrum of a Truncated Sample. For a truncated sample, 
corresponding to a given wy , we define a sample correlation as 


1 fe h h 
(4.1) pr(h, w) = T s(e — fan) f(¢-+4, 0) dé, for |h| 


|h| /2 


WA 


T, 


and 


IV 
_ 


(4.2) pr(h, wo) = 0, for |h| = 
Let us remark that both (4.1) and (4.2) can be replaced by the formula 


oo 
pr(h, wo) = +h fr (8 — fon) fo(@ +5,0) dé 


ie 


=F Srl &, wo) fr(E + | h|, wo) dé, for all h. 


(4.3) 


From H, it follows that the correlation pr(h, w) exists for almost all samples 
(i.e., with probability one). The great advantage of our definition is that the 
correlation pr(h, w) is a positive-definite function of h, uniformly continuous in h (for 
each w for which it exists). 

If we had used as correlation of the truncated sample, as is very often done, the 
function 
prth, w) 


_ lA 
1 TT 


(4.4) prlh, w) = 
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then we would have completely missed these important properties of pr(h, w). 
Indeed, jr(h, w) would have been, in general, discontinuous at |h| = T and no 
longer a positive-definite function. Thus j7(h, #) would not have been the Fourier 
transform of a function ¥7(A, w). 

The only advantage gained by using fr(h, w) is that if f is constant, jr is also a 
constant while pr is not. This fact was probably the reason which led the statisti- 
cians to use this definition for the correlation of a truncated function. However, the 
nonexistence of a spectrum ¥r(A, w) which may have to be reintroduced later by 
various artifices, may lead to serious complications in the estimation of the correla- 
tion functions, particularly when numerical methods are used for that purpose. 

The spectrum ¥r(A, w) is very simply connected to the complex Fourier trans- 
form of the sample 


T +2 
(4.5) ar(A,w) = i [ e f(t, w) dt = 4 | e™ fe( t, w) dt. 
T 0 r 00 
Due to (2.3), the Fourier transform exists for almost all samples and, by Plancherel’s 


theorem, [8, Vol. 2, p. 742], ar(A, w) € L'[—«#, +] (but not, in gen- 
eral, to L[— ©, + ©]). From (4.5) we have: 





+o 
(46) lara) [= a [EP felt, ) fas, «) at ds. 
Let us consider the (¢, s) plane and make the change of variables 
t=§— 4 s=i+ : . 
(4.7) 
t= i ; h=s-—t. 


2 
For any F(t, s) € L(R’) 


[[ res) aas a IE r( —Aie+h) aca 
8) 


CUE r@-$: +4) ae] an 


the last formula being a consequence of Fubini’s theorem [8, Vol. 1, p. 631]. 
Using this transformation, we can write (4.6) in the new form 


(4, 


(4.10) | ar(d, w) |? = ‘ tt e or(h, w) dh. 
Using the fact that pr(—h, w) = pr(h, ») we have 

(4.11) baeGhi, on) fo rs [ sgt wo) ann weil 
Thus we obtain immediately 


(4.12) ites oid -2 [ ith wi) nen ted. @ F | ar(, w) 
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which gives the expression of the spectrum in terms of the Fourier transform of the 
truncated sample. 


Let us note the following properties of this function 


(4.13) ¥r(A, w) Pa 0, 
(4.14) ¥r(—A, w) _ ¥r(A, w), 
(4.15) ¥r(A, w) € L(O, +). 


Due to this last property we can invert the Fourier transform (4.12) and we obtain 
the reciprocal formula 


+20 


(4.16) pr(h,w) = F r(A, w) cos Ah dd. 


5. Sub-Correlation and Sub-Spectrum of a Random Function. Let us now define 
the sub-correlation R;(h) of the random function f;(t, w) as the average of the 
sample-correlations pr(h, w), i.e., 


(5.1) Rr(h) = prlh, w). 
Obviously, we have 
(5.2) Rr(h) = 0, |h| > T. 


and for all other values of h 


* h h 
(5.3) rein) = [[F | §(e-hw)s(e+h,u) ae] du, |h| ST. 
OLT Jiaye 2 2 
Inverting the double integral in accordance with Fubini’s theorem we find 
. rT h h 
(5.4) Rr(h) = = r(z-he+ hax Osh ST. 
T J iayie 2 2 


Let us refer to the change of variables (4.7) and let us consider the (¢, s) plane 
(Figure 1). We introduce the following notation 


(5.5) Ar = {(t,8):0 sts T, 0sss T}. 
Then obviously 


(5.6) Rr(h) = a X Integral of I'(t, s) along the segment AB of 4(h) contained 


in Ar. 


If we define 

r's(t, 8) = I(é, 8), (t,8) € Ar 
(5.7) 

rr(t, 8) vee 0, (t, 8) € Ar 
we can also write (5.4) in the following way 

© xs h h 


which is true for all values of h, giving Rr(h) = Oforh < T orh > T. Obviously, 
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Fig. 1.—Change of variables. 


one has also 


h) =F] rele + [hI ae 
(5.9) 


—_ 


T—|h| 
iy T(z, & + |h|) ak. 


Thus, to compute Rr(h) for all h it is sufficient to know the covariance T(t, s) in the 
square Ar. 

Let us define the sub-spectrum of the random function f(t, w) as the average over 
the spectra of the truncated random function 


(5.10) er(A) = ori, w). 


Then we have the two reciprocal formulas 


+00 
(5.11) Rr(h) = [ ¢r(X) cos Ah dd 
0 


(5.12) eit) a =[ Rr(h) cos dh dh. 
FT 


From (4.13) it follows that g7(A) 2 0; thus, by S. Bochner’s theorem [2] Rr(h) is a 
continuous correlation function. 


6. Correlation Function of a Random Function. We define the correlation function 
R(h) of the random function f(t, w) by 


T—|h|/2 
(6.1) -R(h) = lim Re(h) = lim 5 r(:- h +4) a 


7 g 
toto 1 Jinj/e 3° 


- 


if this limit exists for every real h. 
If the random function is stationary (wide sense), then by definition 


(62) ris) =o(s—0, r(e-h,e+*) = om), 





an 


(6 


TI 


(6 
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and 
Re(h) = (1 J oy) o(h), [hl <7, 
(6.3) T 
Rr(h) = 0, |h| 2 T. 
Thus the limit of Rr(h) exists and 
(6.4) R(h) = lim (1 a Li) p(h) = p(h). 
T-2 y 


Hence, for a stationary random function our definition gives the classical result, 
but we can also apply (6.1) to non-stationary random functions. 
Let us consider as an example the random function 


_ Wit) 
vt 
where W(t) is the classical Wiener-Lévy function, giving the abscissa at time ¢ of a 


particle, starting from the origin at time ¢ = 0, and subjected to one-dimensional 
Brownian motion. This function is certainly not stationary; it has the covariance 


S(t, #) 


(ts) = 4/-, O<t 


IA 
) 


IIA 


asad 
T(t, s) N 7? O<sSt. 


According to our definition this nonstationary random function has the correlation 


ak 
R(h) - heel rf WV 41h] dé = a for all Ah. 


7. Spectrum of a Non-Stationary Random Function. As far as the spectrum is 
concerned, ¢7r(A) does not, in general, tend toward a limit when T — + ©, even if 
the correlation R(h) exists, but exactly as in the stationary case [10, Vol. 2, p. 
164-166] it can be shown, using Paul Lévy’s continuity theorem [12, p. 195], that 
the integrated spectrum 


» 


(7.1) Sx(a) = | er(n) dn 
0 

does in fact tend toward a limit 

(7.2) S(A) = lim S,7(A) 
T-+0 


if the correlation R(h) defined by (6.1) exists and is continuous. (This is not 
necessarily true; R(h) being the limit of a sequence of continuous functions can be 
discontinuous.) Thus, 
(73) R(h) = [ cos xh dS(a), 

0 


this being a Fourier-Stieltjes integral and the spectrum S(\) being a monotonic 
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non-decreasing function such that S(0) = 0, S(+«) < +. In general, S(A) is 
discontinuous and the spectrum has a countable number of lines corresponding to 
a finite amount of energy. 

The spectrum S(\) can be computed from the correlation R(h) by applying 
the Paul Lévy inversion formula [12, p. 166] 





or ° 
(7.4) id ~ din ae f Se a) aw, 

A>+o T h 
which is valid if \; and ), are continuity points of S(A). When R(h) is known the 
spectrum is thus defined for every \ > 0 with the exception of, at most, a countable 
number of discontinuity points. 


8. Estimation of Correlations for a Non-Stationary Function. For a non-stationary 
random function, even if we have not only truncated samples of the function, but 
also its covariance in A; , this does not give us sufficient information to determine 
R(h). It is obvious, from (6.1), that large values of are most important in deter- 
mining R(h) (even at small values of h). The knowledge of I(t, s) in the square 
Ar only does not give us any information about its values for large — on 5(h) (See 
Fig. 1). 

We shall consider here one class of random functions which is not stationary, 
but on which information is given, which enables us to make an estimate of R(h). 

This class is defined by the condition that T (« - 4 § + *) is periodic in &. This 
condition means that the covariance is invariant under a translation 7 parallel to 
the first bissectrix (Fig. 1) 


(8.1) r(i¢+7,8+ 7) = P(t,s). 


The scope of the implications of this hypothesis for applications is better under- 
stood if one points out that (8.1) is satisfied when the phenomenon represented by 
the function f(t, w) is the result of the superposition of two phenomena, one sta- 
tionary and the other periodic. Thus 


(8.2) f(t, o) = filt, o) + falt, w) 
where f,(t, w) is a stationary (wide sense) random function: 
(8.3) Ait, w)fis, ©) = a(t — 8) 
and f2(t, w) is periodic 
(8.4) felt + 7,0) = folt, w) 
fils, w)falt, o) = T2(t, 8) 
I(t + r,s) = T(t, + 7) = P(t, 8), 


(8.5) 


the two random functions f, and f. being uncorrelated 


(8.6) filt, w)fo(s, w) = 0. 





Th 
(8.7 


and 


(8.8 
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Thus one has 
(8.7) rt, 8) = pilt oe 8) + l(t, 8) 


and I(t, s) satisfies (8.1). 
In this case the limit (6.1) obviously exists and the correlation is given by 


(88) R(h) == [ rg e+ [hI de 

The sub-correlation for T = N7r(where N is an integer), due to (5.9), is equal to 
1 Nr—|h| 

(8.9) Ry(h) =~ [TE + (hI) ae 


As a consequence of the periodicity 


‘ 1 Jal 
R(h) — Ruch) = wil r(é,é + |h|) dé. 


We thus have the upper bound 
1 pil 
(8.10) | Rh) — Reh) S-[ IT E+ IAD Lae 
T 


for the error obtained by using the sub-correlation Ry,(h) instead of the correlation 
R(h). The upper bound (8.10) is a function of h; however, at a given h, this bound 
tends to 0 as 1/N. 

Let us consider the case when (k — 1)r & | h| S kr (k integer). Using the pe- 
riodicity of and Schwarz inequality we have 

38 4 

(8.11) | R(h) — Rwe(h) | $= [ TUE 8) ae 
When the covariance I(t, s) is known for one period over the diagonal t = s we can 
compute 


1 Tv 
(8.12) A=-[ r(,8) de. 
T J 
Finally, the upper bound for the approximation is given by 
(8.13) | R(h) — Rw-(h) | < AS, (k —1)r Sh Ske. 


Let us now consider Rr(h) when Nr < T < (N + 1)r. We have 


ll 


. T—|h| 
Rs(h) — FF Rel) = Gf MRE + LAD a 


1 T 
pL PEt IAL a 


Thus, due to the periodicity of I, 


(8.14) | Roth) — X7 R,.(h) | 


< 
T 4 


<4 
N° 


| 
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From (8.13) and (8.14) we obtain finally 





Rh) — Nt R(h) sai t* Nec T<(N +1), 


(8.15) ie 


(k —1)r Sh < kr. 
This relation shows that the approximation is good if k/N is small. 


9. Accuracy of Estimates of Correlations. In estimating the correlation R(h) of 
a random function from truncated samples there are two steps: 

(a) From the correlations pr(h, w;) of the truncated samples we estimate the 
sub-correlation R7(h), 

(b) From R7(h) we compute R(h). 

In the preceding section for the case of a random function with a periodic covari- 
ance we have solved problem (b) at large values of 7. Now, let us look at the prob- 
lem (a), namely, how to evaluate the approximation with which one determines 
R,(h) from the average of the correlations p7(h, w;) of a number q of samples 

. 1 = 
(9.1) Rr g(h) = ; > prth, w;). 
= 
It appears that the best way to make such an evaluation is to determine the vari- 
ance of the random variables pr(h, w;) about their mean value R,(h); if this vari- 
ance is small enough we can expect that for a reasonably large number g of samples 
the estimate Rr(h) will be fairly good. 

In order to compute this variance, in addition to Hypotheses H, to Hy, of Section 
2, we shall assume that 

H; . f(t, #) € L*(Q) for all ?’s. 

This insures the existence of the fourth-order moment 


(9.2) Ot, tr, te, ts) = [ flty, w)f (te, o)f(te, w)f (ty, o) dy 


for all [t; ,4,h, tJ. 

We shall, moreover, assume that 

Hy, . W(t, b, ts, ts) € L(A) for every finite parallelepiped oi the four-dimen- 
sional space (t; , te , ts , ta). 

Let us observe that the fourth-order moment exists in the important case of 
normal random functions, i.e., when the n random variables, f(t, , w), --- f(t. , ) 
[t: , ---+ , tr] arbitrary follow an n-variate normal (Gaussian) law. 

From (4.3) and (5.8) we find that the departure of a correlation for a particular 
truncated sample from the mean value taken over the correlations for all samples is 
given by 


oo 
pr(h,w) — Rr(h) = of [re(: —he+h) 


— fe (¢- 2,0) se(¢+8,0) | ae 





Aft 
(9. 


wit 


In 


(9 


th 


of 
») 


ar 


de. 
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After taking a square and averaging, we obtain 


(93) ox(h)* = [prliga) — ReH)P = 7 [[ Fe(e,0,h) dé dn 





with 


h 
F7(é, n, h) > mr(—hethin—tin th) 


h h h h 
—re(e—Be+h)re(n— sat). 


In the above equation I'7(¢, s) is defined by (5.7) and 9M, by the two relations 
Mrli,&,&,4) = M(h,b,b, &) 
when (4,4) € Ar and(t,,%) € Ar 


(9.4) 


(9.5) 


and 
(9.6) Mr( ty : te : tz 9 ts) = 0 when (4; : t) é Ar or (ts ’ ts) € Ar. 


Equation (9.4) shows that whenever the covariance T and the fourth-order moment 
9 are known, we can evaluate the variance, or(h)*, of the random variables 
pr(h, w) about their mean value R,(h). In particular, for h = 0, we have 


(97) x(0)* = 75 [ff taw(g, & ny.) — PCE, OT Ca, 9] a a 


10. Random Step- Function with Periodic Covariance. As an example let us take 
the random function defined as follows* 


(10.1) V(t,w) = X,(w) n-—1lst<n,n = 1,2,--- 
where 
(10.2) Xi(w), +--+ , Xn(w), --- 


is a sequence of independent random variables, taking only the values —1 and +1 
with equal probability 


(10.3) Prob [X, = —1] = Prob [X, = +1] = }. 


The random function V(t, w) is essentially the same as a function considered by 
Norbert Wiener in his pioneering work on correlation and spectrum, [15, p. 151] 
except for the fact that his function is defined for —« < t < + and ours only 
in [(0, +]. 

From (10.3) it follows that 


(10.4) XxX, =0 
(10.5) XaXn = dan, 

* As measure space Q, one can take the interval [0, 1],0 < w < 1, with X,(w) = 2a, — 1, 
where a, is the nth coefficient in the binary development w = a,/2 + --- + a,/2" + --- . The 


measure u on [0, 1] is the Lebesgue measure p{w: a, = 1} = u{w:a, = 0} = f. 
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where 6,,,, is the Kronecker symbol. If one defines 
Amn» = {(t,8):m—-lst<mn-—-1358 <n} 
then 


n=+o 


M(t,s)=1, (t,8)€ U Aw, 


n=l 


n=+o 
T(t,s)=0, (t,8)€ UA. 
n=l 
Obviously, the random function V(t, w) is not stationary, but its covariance, when 
put in the form T (« oe =. e+ *) is periodic, with period 1 in ¢ (see Fig. 2). Thus 


the results of Section 8 apply, and the non-stationary random function V(t, w) has, 
according to (8.8), a correlation 


R(h) = (1 — |h)), |A| £1, 
R(h) = 0, |h| 21. 
Let us first consider the case T = N, where N is an integer. We find that 

pw(h, w) = (k — |h|)¥w(k — 1,0) + (1 — k + |h|)¥x(k, @), 
(10.7) kK-—-1s|h| Sk, l1sksWN —-1, 
py(h, w) = 0, JA, =N, 


(10.6) 


IV 


where Yy(0, w), Yw(1, w) --- Yw(N — 1, w) represent the random variables 


j=N—k 


(10.8) Yu(k,0) =< DL XAw)Xsa(w), k = 0,1, --- (NW — 1) 
N jal 


4 4S 
S) wet 


3 ho 
me yy 
5 bt a Y) 

















Ny Ny 
0 | 2 5 


Fic. 2.—The (¢, s) diagram for a random step-function. '(é — ‘, -+ r) is periodic in ~ 
and for a translation parallel to ¢ = s, (t+ 1,s+1) =T(t,s). 
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in particular 
(10.9) Yy(0,) = 1. 


For each sample w, the sample-correlation py(h, w) is represented by a polygonal 
line (Fig. 3). The ordinates Yy(1, w), Yw(2, w) --- Yw(N — 1, w) corresponding 
toh = 1, 2, --- (N — 1), are the random variables defined by (10.8). We see 
immediately that Yy(k, w) can only take the values 





_n=s N-—-k-2 N-k-—2 N —k 


I ote i) mel eile 





following the binomial law. As a consequence 


(10.10) | ov(k, a) | S1—<, Tet een Pe 


We obtain for the mean, variances, and covariances of the ordinates Yy(k, w) 
respectively, 
(10.11) Yx(k, ow) = 0 k =1,2,---,(N — 1) 


N-—k 


(10.12) Yy(k, w)? = — Tr? 







ye RU) = Ry 
¥, (rl \\ 


Pylhw) 


a. Sw 
\ Y ba \ ie ee 
uy * i \- 
1\ 1 2 3 N-I N h 
Y (1) 


a 


Fig. 3.—The polygonal line pw(h, w) represents a sample-correlation of the step-function. 


Ry(h) = pw(h, w) is the sub-correlation, which for the step-function, is equal to the correla- 
tion function R(h). 
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and 
(10.13) Yx(k, w)Yyx(l, wo) = 0, k <1. 
From (10.7), (10.9) and (10.11) we deduce 
(A, =1-|h ? h <= 1 

(10.14) pw(h, w) | A | | A | 

pw(h, w) = 0, |hi21 
Thus 
(10.15) Ry(h) = pwlh, w) = R(h). 

Let us next suppose N < T < N + 1. We find 
N k—-N+h T—k—-h 

(10.16) pr(h, w) is 4 T pw(h, w) + ea Xi X was + 7 Xi Xy ’ 


N-k<h<T-—k, k= 1,2,---N. 
N T—wN 
(10.17) pr(h, w) ca T py(h, w) sr Xy X wit; 
T-—-k<h<N+1-k, k=1,2,--- N. 
ThusifN <T<N+1, 


(10.18) —_|pr(h,w) — py(h,w)| (1 is *) low(h, «)| + ( ‘a *) 


As a consequence of (10.10), we have 


2 
Pee -- y 
(10.19) lpr(h, w) — py(h,w)| < N41 for NST<N+1. 


For large values of N the right-hand side of this inequality is as small as we want. 
Thus we can always, in studying pr(h, w), suppose that T has an integer value N. 

Let us now compute the variance oy(h)’ of the sample-correlations about their 
mean value. We will make the computation only for 7 = N. Due to (10.14) one has 


on(h)” = py(h, w)?. 
From (10.7), (10.12) and (10.13) we obtain 


ah +2 
N? 


+(1—k+[ajpr—F, b+ 1s |b SE, Bm 1,2,---N 1. 


pu(h, w)® = (k — |h|)° 





(10.20) 





In particular for h = k 
(N — k) 
ies a 
The proposition that the random variables py(h, w) tend toward their limit R(h) 
for almost all samples, 


Prob [lim py(h,w) = R(h)| = 


N>+0 


(10.21) on(k)” = py(k, w)? = 





wil 


(1 


on 


a Th 


es lm Oe OemlCO DOCH 


nt. 


eir 
as 


h) 





CORRELATIONS AND SPECTRA FOR NON-STATIONARY RANDOM FUNCTIONS 15 


will now be established. First from (10.8) we compute 





(10.21) Yk, a) = (N — = 2k + 1) 


Let us now observe that for a fixed N, when k takes the values 1, 2, --- (N — 1) 
one has Y x(k, w)* < = . Thus, k now being fixed, we have 


N? 
N=+o N=+o0 1 
> Yxrkot S3 > = <+o 
N=1 ya N? 


and, therefore 
N=+0 
Prob | >, Yalkw)* < + = | =1, 
N=1 


Here, we use the following criterion for the almost sure absolute convergence for a 
series of random variables. If }>{* | X, | < +, then Prob [}-{*| X, | < + @] 
= 1. This criterion applies even when the random variables X, , --- , X2, --- are 
not independent. However, the convergence of the series implies that Yy(k, #) 0. 
We have thus proved that 


Prob [| lim Yy(k,#) = 0] = 1 
N>+0 


for each fixed k = 1. 

Let us now take a fixed interval |h| <= M, where M is an integer. In thisinterval, 
there are exactly M points, h = 1, h = 2, --- h = M, which completely determine 
the polygonal line. Each of the M ordinates Yy(1, w), Yw(2, w) --- Yr(M, @) 
tends toward 0 with probability one. Their number being finite, this evidently 
implies 

Prob [| lim Yy(k,w) = 0 for k = 1,2,--- M|=1. 


N7>+2 


We have thus proved that 
(10.22) Prob [ lim py(h, #) = R(h)] = 1 
N+ 


in the finite interval | h | < M and the proof is complete because M could be taken 
arbitrarily large. 


11. Correlation Estimates for a Continuous Step-Function Constructed Using 
Sequences of Random Numbers. A continuous step-function has been constructed 
using a table of random numbers which are listed in 100 groups of 1000 digits each 
[11]. For our analysis we selected the 200 first digits of each of the 100 groups and 
divided them into five consecutive segments of 40 digits. All the even digits were 
then replaced by +1 and the odd digits by —1. We thus obtain one set of 500 
sequences of 40 digits (+1, or —1). Other sets were constituted by combining two 
or more consecutive segments of 40 digits. As a result we have obtained the follow- 
ing five sets of experimental functions without overlapping sequences within each 
set: 500 sequences of 40 digits; 200 sequences of 80 digits; 100 sequences of 120, 
160 and 200 digits, respectively. 
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Let us consider a set of q sequences with N digits in each sequence. A sample- 
correlation for the sequence is determined for T = N by (10.8). The average taken 


over the q sample-correlations of the set determines the estimate for the sub- 
correlation, cf. equation (9.1). 


(11.1) Ry (k) = - x, sith. ox), 


In a sequence of N digits there are N — k products X ,(w)Xj4.(w) equal to 
either +1 or —1. If ny(k, w) is the number of products equal to +1 then 


(11.2) pu(k, w) = = (N — k — 2ne(k, o)] 
and 
(11.3) Ry,(k) = 2 


q 
= —_ Na > ny(k, wi). 


As it has been shown in Section 10 a sample-correlation py(h, w) is represented 
by a polygonal line. The vertices of this polygonal line correspond to h = k = 
0, 1, 2, --- , (N — 1) and are given by (11.2). Similarly, the estimate Ry,,(k), 
for the sub-correlation is represented by a polygonal line which is determined by 
(11.3). It is, therefore, sufficient to determine the values of py(k, w) and Ry,.(k) 
at the vertices of the polygonal lines to have the corresponding sample-correlations 
pw(h, w) and the estimate for the sub-correlation Ry,,(h) for the continuous step- 
function. Figure 4 illustrates several examples of sample-correlation py(h, w) for 
individual sequences of digits. Numerical data for Ry,.(k), for the five sets of 
sequences, are listed in Table I and a few of them are illustrated on Figure 4. The 
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Fic. 4.—Experimental examples of sample-correlations py(k, w) obtained for a random 
sequence of digits —1, +1, are represented using light interrupted lines. The heavy polygonal 


line represents the estimate of the sub-correlation Ry,,(k) obtained by averaging over g = 500 
sample-correlations. 
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theoretical sub-correlation Ry(h) is in the present case equal to the correlation 
R(h) (see (10.15)) and, therefore, the experimental results for Ry,.(k) should be 
compared to the theoretical value 


(11.4) Ry(k) = R(k) = 0, k21. 


The departure of the value for a sample-correlation from the estimated sub- 
correlation, as determined from g samples, is given by 


pw(k, w) — Ryo(k). 
The variance oy,.(k)” for such departures is given by the relation 
ona(k)? = ~¥ lowly) — Ry oH 
(11.5) wis us ‘ 
- We ) » [ -na(t, w) + a ny(k, ») | ; 
Numerical results for the experimental variance oy,.(k)’ are listed in Table II and 


compared with the theoretical values, ox(k)*, computed by using the relation 
(10.21). 


TaBLeE I 
Estimates for the Sub-Correlations of Sequences of Random Numbers 














Yee tee 40 | 680 | 120 160 | 200 
” beahd heh ohh Sek. 500 | «62000 | «=100 | 100 | 100 
Rx..(1) | —0.0041 | —0.0026 | —0.0047 | —0.0023 | —0.0053 
Ry¢(2) | +0.0053 | +0.0040 | +0.0143 | +0.0033 | +0.0054 
Rx..¢(3) | +0.0066 | +0.0071 | +0.0087 | +0.0081 | +0.0073 
Ry (4) | +0.0083 | +0.0036 | +0.0095 | +0.0061 | +0.0072 
Taste II 


Comparison between Experimental Variances ox ..(k) and Theoretical 
Variances ox(k) 




















_ Pe ae ee | 40 | 8 | 120 | 160 200 
TAI BE PA Le 3. | 500 | 200 | 100 | 100 100 
| | | 
ow .¢(1)? | 0.0253 | 0.0125 | 0.0094 | 0.0074 | 0.0066 
on(1)? | 0.0244 | 0.0123 | 0.0083 | 0.0062 | 0.0050 
ox o(2)? | 0.0231 | 0.0125 | 0.0076 | 0.0061 | 0.0048 
on(2)? | 0.0238 | 0.0122 | 0.0082 | 0.0062 | 0.0050 
ox. (3)? | 0.0232 | 0.0121 | 0.0074 | 0.0059 | 0.0043 
ox (3)? | 0.0231 | 0.0120 | 0.0081 | 0.0061 | 0.0049 
ov .«(4)* | 0.0198 | 0.0089 | 0.0065 | 0.0038 | 0.0035 
ox (4)? | 0.0225 | 0.0119 | 0.0081 | 0.0061 | 0.0049 
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15. Frequency Distributions of Sample Correlations. Let us denote as 
q5 lnw(k, w) 4 ul] 


the number of sequences of a set for which ny(k, w) is equal to u. We find then, for 
the frequency distribution of sample-correlations 


(12.1) Falew(k, w) = U] = F[nw(k, w) = ul, 


where 
U = 5 (N —k—2u). 


Numerical values for the frequency distribution of sample-correlations can easily 
be obtained from Table III where gS,[ny(k, w) = wu] is listed for the five sets of 
sequences. It may be noted here that (11.3) can also be written as 


: N-k 
(122) i ays 2. a > eticfn(ho) = al, 


u=0 


which is often more convenient for numerical computations. 












































TaBLe III 
Frequency Distributions of Sample-Correlations Fq(py = U) = Fq(nw = u) 

500 Fs00 (n4 = u) 200 $200 (nso = u) | | 100 F100 (ni = u) 
nwo 180 | | mi | : --. 
k=1k=2;)k=3\k=4 k=l\k=2k=3|k=4| k= 1k =2\k =3)k=4 
—_—|-_——— —_—_—_—|————— ———} —— | —— — | | ‘al fees ee tee i-——— 
o| of of of 1] |a} of of 1] o| |) 0}1] 0] o 
10} 0} O} O| 1 %| 0} O| O 1] /4/1/0|]24| 2 
11] 4] 3| 2 3 2%/ 0} O| O 0| |47;0}1]21] 0O 
_s wiie) s | 8 27; O| 1 1 | 0] |48;}1]|2] 1 1 
13} 8| 16] 10 | 19 23| 0]; 1} 2] 0| |49}2|2)]1 0 
14) 14| 13] 18 | 23 2; 2} 1] 1] O| |50/ 2] 0] 3 1 
15 | 23 27 | 49 0}; 0; 3] 3] 1] |S} 1)]3);04] 2 
16| 30| 45| 26 | 44 31) 4) 3) 2) 3 | | 52) 3 oe ee 
17| 49| 53) 46 | 61 32; 4/ 8| 7] 8 53| 6 | 5 | 4 7 
18| 50} 51| 55 | 74 33} 9| 7| 12 | u | js] 2] 3 | 4 7 
19| 58| 66| 69 | 72 34/ 5| 12] 14 | 13 | |55| 7 | 7 | 10 9 
20| 56| 70| 68 | 65 35| 13| 14 2 | 19 | |6| 5 [13 | 8 8 
21| 68| 57| 59 | 38 36| 18] 15| 13 | 18 | |57) 6/7) 9] 11 
22| 51/ 20) 46 | 18 37 | 11) 13) 15 | 18 | | 58) 7 | 7 | 9 | 10 
23 | 38| 32] 2 | 6 38} 16| 15] 19 |. 18 59 | 5 = ee. 
24| 22} 9| 20/| 8 39| 16| 20| 20 | 29 | |6|7/6|8]| 4 
25| 11| 12] 13 | 5 | | 40] 18] 10 20 | 23 | |61/ 6 | 8 | 2 4 
%| 7| 4/ 8 | 0 41/ 14] 11} 15 | 6 | |62} 5 |71] 5 3 
27; 2] 3| O 1 42| 19] 21| 12 | 14 | | 63} 4/1 | 6 4 
| 6| 1} 3| 0| | 4 16} 12} 8 | 5 | |64| 9 | 4] 4 4 
| 144; 6] 11) 5 | 5 | |65/ 6 | 2] 3 2 
45; 8| 12; 7 | 0 | || 6/3]51] O 
| 46/ 6| 4] a 4] |a| } 3 | 0] 0 
47; 8| 3| 3 1 | |68/ 0 | 2] 0 0 
48) 3| 2) 1 | 1 | | 6 3/1]0 2 
49) 2) 0/ O | O} |7%| 2) 0/0) 2 
| |}50} 1) O| O} 2] |7}o0/o0}0)] 1 
| 51] O| 1| 1] 0 | |v} 2] 0/11] 0 
| | |52{ O| O| O} O | | 7%] 0] 0 | O 0 
|} |53; O; O} O| O | | 74) 0 a eS es. 
| eS co: ies on ee ae ee ee 

| | | Awe | | 
| | |77/ 1} 0]0| oO 

| 








| 
| 
| 
| 
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TABLE III (Continued) 


100 Fi00(n160 = u) 
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For large values of (N — k), this binomial distribution can be approximated by 


the Gaussian distribution 


relation 
(12.3) 


where 





1 


(N —k — 2u) 


(12.4) 





or 


a) | rrp | 


u = 0,1,- 


2(N — k) 
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Fig. 5.—Experimental frequency distributions of sample-correlations, F,fon(k, w#) = U} 
(full line)are compared with the theoretical binomial frequency distribution (dashed lines). 


In particular, the probability of py(k, w) being equal to the sub-correlation 
Ry(k) = py(k, w) = 0, k 


(which, in the present case, is also equal to the correlation R(k)) is given by 


IV 


, 1 N —-k : 
(12.5) Prob low (k, w) - 0) i Qn-k (it 2. ib) ’ k21 
with the approximated formula 
2 
(12.6) Prob [py(k, w) = 0) = ———————. , k>1 
’ /2x(N — k) 


In Figure 5 we are comparing some of the experimental frequency distributions 
with these theoretical results. 
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On Numerical Integration of Ordinary 
Differential Equations 


By Arnold Nordsieck 


Abstract. A reliable efficient general-purpose method for automatic digital com- 
puter integration of systems of ordinary differential equations is described. The 
method operates with the current values of the higher derivatives of a polynomial 
approximating the solution. It is thoroughly stable under all circumstances, in- 
corporates automatic starting and automatic choice and revision of elementary 
interval size, approximately minimizes the amount of computation for a specified 
accuracy of solution, and applies to any system of differential equations with deriva- 
tives continuous or piecewise continuous with finite jumps. ILLIAC library sub- 
routine #F7, University of Illinois Digital Computer Laboratory, is a digital 
computer program applying this method. 

1. Introduction. A typical common scientific application of automatic digital 
computers is the integration of systems of ordinary differential equations. The 
author has developed a general-purpose method for doing this and explains the 
method here. While it is primarily designed to optimize the efficiency of large-scale 
calculations on automatic computers, its essential procedures also lend themselves 
well to hand computation. The method has the following characteristics, all of 
which are requisite to a satisfactory general-purpose method: 

a. Thorough stability with a large margin of safety under all circumstances. 
(Instabilities in the subject differential equations themselves are, of course, re- 
flected in the solution, but no further instabilities are introduced by the numerical 
procedures.) 

b. Any integration is started with only the essential initial conditions, i.e. 
there is a built-in automatic starting procedure. 

c. An optimum elementary interval size is automatically chosen, and the choice 
is automatically revised either upward or downward in the course of an integration, 
to provide the specified accuracy of solution in the minimum number of elementary 
steps. 

d. The derivatives need be computed just twice per elementary step, which is 
the minimum consistent with controlling accuracy. 

e. Any system of equations 

OU = fila, ty Ys -*) 62 1,2,---n 
(1) 
dy 
(often written a f(a, y) for short) 


can be treated for which the f; are either continuous or piecewise continuous func- 
tions with finite jumps. 

f. The solution is computed at (although not necessarily only at) equally spaced 
values of the independent variable x, with specifiable spacing. 


Received May 6, 1961. The research presented in this paper was supported in part by the 
U.S. Army, Navy and Air Force. 
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Further useful though perhaps not indispensable characteristics of the method 
are: 

g. Enough numerical information is developed to make interpolation or evalua- 
tion of functions (e.g., roots) of the solution possible with accuracy equivalent to 
the solution accuracy. 

h. The sense of integration can be reversed. 

Characteristic a) is essential for getting trustworthy results in lengthy auto- 
matic computations because the number of elementary steps may be as large as 10° 
or 10° or more, and disturbances in unstable methods typically grow exponentially 
with the number of steps. Characteristic b) is not only a convenience but also 
insures that in the integration of intrinsically unstable equations, in which early 
errors tend to be strongly magnified, the starting errors do not dominate. Charac- 
teristic c) relieves the human being of the often difficult task of determining the 
correct interval in advance. Where the human being must specify the interval for a 
computation not to be performed by himself he tends to make up for uncertainty by 
a conservatively small interval choice. Characteristics c) and d) together thus make 
for efficient use of computer time, and the saving in computer time can easily be a 
factor of 10 or even much more in the handling of problems in which the interval 
should vary. 

In regard to the question of relating our method to previously available methods, 
we wish to make clear at the outset that it is equivalent to a reformulation of the 
method of Adams [1, p. 53-55], [2, p. 81-82] for it uses effectively the same quadra- 
ture formula as does Adams. However, the formulation and the point of view are so 
different that it is instructive and seems appropriate to explain the method starting 
from first principles, as we shall do below, rather than starting from Adams’ quad- 
rature formula. 

Presently available methods may be divided into two classes: those involving no 
memory and those involving some memory, of the past behavior of the solution. 
The Runge-Kutta methods [1, p. 72-75], [2, p. 59-74] are typical of the first class, 
the Milne methods [1, p. 64-70], [2, p. 84] and the Adams methods of the second. 
It has been clear for some time that the methods with memory are superior in 
accuracy for a given elementary interval size and a given amount of computational 
labor since they permit a better approximating curve to be fitted over the elemen- 
tary interval. Our method involves such memory. In return for this superiority of 
the methods with memory we must cope with two problems quite foreign to the 
memoryless methods: how to start off, since at the beginning there is nothing to 
remember; and how to prevent the remembered numerical information from behav- 
ing unstably. 

Two further problems must be dealt with in order to implement the automatic 
choice and revision of the elementary interval, namely, choosing which quantities to 
remember in such a way that the interval may be changed rapidly and conveniently 
and developing an appropriate set of rules for controlling the interval size. Thus the 
four major problems are: automatic starting, stability, choice of quantities to 
remember and interval control logic. The last of these four is the most intricate. 

As with most methods, there exist lower and higher “order” versions of this 
method. The author prefers to use the term ‘‘degree” rather than “order”, since all 
methods are ultimately equivalent to finding a polynomial of some given degree 
approximating the solution of the system of equations, and since the term “order” 
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is already standardized usage for the number of equations n in (1). We have chosen 
and recommend degree 5, which corresponds to a truncation error 0(h’) per elemen- 
tary step of length h, for large-scale digital computer operations. This represents an 
advantageous return in accuracy per step with quite large steps, while still not 
overdoing the accuracy when the choice of h is limited to inverse powers of 2, as is 
natural in a binary computer. 

The order n of the system (1) is immaterial to a large part of our discussion, so 
that we can advantageously use the simpler notation dy/dx = f(x, y) for (1), 
regarding y and f as vector-like objects with n real numbers as components. The 
independent variable x is, of course, a single real number. Whenever the multi- 
component character of y and f makes a significant difference in the discussion we 
shall so note. 

In Section 2 the choice of quantities to be remembered is discussed, in 3 the 
numerical procedure and the associated stability theory are developed, in 4 certain 
parameters of the method are adjusted for optimum stability and accuracy, in 5 
the procedure for modifying the interval is given, in 6 the characteristic behavior of 
the remembered quantities is described, in 7 error estimation is discussed, in 8 the 
automatic interval control logic is developed, in 9 automatic starting is described 
and finally in Section 10 the results of certain test problems done by this method 
are exhibited. In Appendix A are collected the working formulas and error estimates 
for degrees 3 through 6 of the approximating polynomial. Appendix B contains a 
schematic flow chart for programming the method for a digital computer, with com- 
puting time estimates. Appendix C is a discussion of control of roundoff errors in 
iterative numerical procedures. 


2. Choice of Quantities to Remember. It is immediately clear that quantities like 
differences y(x) — y(x — h), y(x — h) — y(x — 2h), ete., and/or higher differences 
would constitute a poor choice to remember, for changing the interval in terms of 
these is a cumbersome process involving much interpolation and/or extrapolation. 
(Ignoring the remembered quantities whenever the interval is to be changed and 
starting again “from scratch” would entail serious loss of accuracy and of time). 

We take our cue from the remark above to the effect that all methods of numeri- 
cal integration are equivalent to finding an approximating polynomial for 
y(x). Of the many ways of specifying a polynomial of degree m by m + 1 constants 
there is one way which is interval-independent, namely: to specify the Oth to mth 
derivatives of the polynomial evaluated at the current value of x. These particular 
m + 1 quantities specify the same polynomial no matter what the interval is, 
being in fact defined with no reference to an interval at all. They would be ideal 
from the point of view of interval modification. However, they are not suitable for 
automatic computation because the higher derivatives may vary enormously in 
magnitude and are thus not conveniently stored in a “fixed-point” arithmetic 
operation.* 








* The discussion in the present paper is limited to ‘‘fixed-point’’ arithmetic procedures. 
The question whether a ‘‘floating-point’’ version of the method could be made safe against loss 
or illusory gain of significance of the quantities in the course of a long computation, and other- 
wise trustworthy, is for future investigation. The possible freedom to store just the higher 
derivatives of the approximating polynomial and the increased freedom from scaling problems 
certainly suggest that one investigate the floating-point possibility. 
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In order to see how to modify our choice so as to cure the latter diffieulty, we 
consider how the m + 1 derivatives would actually be used in the computation. A 
typical important use is, in the first phase of the integration step from z to x + h, 
to “predict”’ a trial value of y(x + h) from the formula: 


y"(x +h) = y(x) + h{ se, y(x)) + a Ps” (x) 
we h° P. ga | h° P. pore h‘ P. sever 
+ ahs t)+7 Ps + els (x) 
where m has been made 5 and P;(x) = y(z), Ps'(x) = f(z, y(z)), Ps” --- Ps” 
are the 6 aforementioned derivatives of the approximating polynomial evaluated 
at x. Formula (2) is written in the special way shown, with one factor h external to 
the { }, because we may expect f to be computed to full register accuracy on occa- 
sion, which suggests that the remaining terms in the { | be kept to the same 
accuracy; and because for the case of small h and many steps (many successive 
applications of formulas like (2)) we can minimize the accumulation of roundoff 
errors in y by keeping log (| h |’) more places in h{ } than we keep in the{ | 
itself. Formula (2) in the form written then suggests that the appropriate quantities 
to store in the computer registers are, besides the always necessary y(z) and 
f(x, y(x)), the four quantities 


2 
a(zr) = 4 P(x) = b(z) = h pz) 


2! ! 

(3) hr h* 
c(z) - ho P,!""(z) d(x) om — eee 

4! 5! 
We may reasonably expect these quantities to stay within register capacity since an 
appropriate choice of A will just cause the successive terms in the | }{ to decrease 
in magnitude no matter how large the Ps themselves become. Although the 
quantities (3) are not completely interval-independent, they depend on the interval 
in such a simple way that interval change involves merely multiplying each by a 
constant, and in the important practical case of a binary computer and intervals 
restricted to inverse powers of 2 the change is achieved simply by shifting the 
numbers. Formula (3) seems accordingly to be essentially the unique sensible 
choice, at least for a fixed point arithmetic procedure. 

We emphasize that the quantities y, f, a, b, c, d as they exist in the computer 
registers and appear in our discussion are formally defined from successive deriva- 
tives of an approximating polynomial, so that they always exist since an approxi- 
mating polynomial always exists, whether or not the exact solution of the original 
problem (1) has five derivatives. If the original problem involves a discontinuous f, 
the quantities a --- d tend to get large because of that, but concurrently tend to 
get small because of interval decrease, with the overall result that they stay within 
register capacity. While the existence of an approximating polynomial is assured, 
its quality as an approximation of the exact solution of (1) depends on how it is 
developed; in subsequent sections we discuss how to develop it in an optimum way. 


3. Taylor’s Theorem Procedure Modified for Stability. In order to have a com- 
pletely defined integration procedure we must have rules for determining all of the 
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quantities y(x + h), f(x + h), a(x + h)---d(x + h) when given y(z), 
f(x), a(x) --- d(x) and the differential equation (1). (The starting problem, namely 
to determine y, f, a, b, c, d at x + h given only y(x) and f(x) and the differential 
equation, is discussed below in Section 9). Consider first the ordinary Taylor’s 
series formulas terminated at h°, which in terms of a, b ---- read: 


y(z + h) = y(x) + Aif(z) + a(x) + bie) + cz) + d(x) + = e(z)} 


i@+h)= f(z) + 2a(z) + 3b(z) + 4c(x) + 5d(x) + 6e(z) 
(4) a(z+h) = a(x) + 3b(x) + 6c(x) + 10d(x) + 15e(z) 
biz +h) = b(xz) + 4c(x) + 10d(x) + 2e(z) 
c(z +h) = e(z) + 5d(xz) + 15e(z) 
d(a+h) = d(x) + 6e(z) 
Here we have introduced one more quantity e(z) analogous to a --- d, which we 


eliminate forthwith by using the differential equation. The system (4) as it stands 
is incomplete, having one less equation than it involves quantities. But by identify- 


ing the second formula of (4) with f(z + h, y(x + h)) calculated from the differen- 
tial equation, we can eliminate e(x) and get: 


y(z + h) = y(z) + Alf(z) + a(z)+ bz) + cz) + d(x) 
+ &Uf@+h, y+ h)) — fr} 

SJ(za+h) = J(z) + 2a(z) + 3b(z) + 4e(x) + 5d(z) 
+ 1Uf@t+h, y@+h)) — fr) 

a(jz +h) = a(x) + 3b(x) + 6c(x) + 10d(z) 
+P U@ + h, y@ + h)) — fr) 

6) b(2 +h) = b(x) + 4c(x) + 10d(z) 
+ *P [f@ + h, ye + h)) — fr) 

e(z +h) = c(z) + 5d(z) 
+ ¥ (f@ +h, ye + h)) — fr) 

d(z +h) = d(x) 


+ 1[f(@t+h, y@t+h)) — fr) 


where f? = f(x) + 2a(x) + 3b(x) + 4c(x) + 5d(x), the “predicted” value of 
f(z +h). 

Now the system (5) augmented by the differential equation is complete, for the 
first equation of (5) and the differential equation together constitute an implicit 
system determining y(z + h) and f(x + h); the second equation of (5) is an identity 
and the next four then determine a(x + h) --- d(x + h) straightforwardly. 

Having arrived at the scheme (5) quite directly from Taylor’s theorem we enter- 
tain the possibility of using it for numerical integration. A small amount of hand 
computation using (5) establishes that it is: a) very accurate indeed, and b) very 
unstable indeed, with small disturbances growing approximately as (—10)° in s 
steps. 
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These two phenomena are closely related. The high accuracy derives from basing 
the scheme directly and exactly on Taylor’s theorem; however, just because it is so 
based it has another property, namely reversibility. If we apply (5) to go from z to 
xz + hand reapply (5) with reversed h to retrace from x + h to z, we recover the 
original quantities y, f -- - d precisely. Now a process reversible in this sense cannot 
be stable, for it cannot damp out small disturbances (i.e., “forget’’ or “lose informa- 
tion’) as it must to be stable. Stated in terms of the eigenvalues of the stability 
matrix M discussed later, reversibility implies that the matrix for backward integra- 
tion is the inverse of the matrix for forward integration, which is inconsistent with 
the condition for stability, namely that for both these matrices all eigenvalues 
except one must lie inside the unit circle. (The only exception to the last statement 
occurs when the stability matrix is 1 X 1, which corresponds to the trapezoidal 
method m = 1 with no “memory.”’) 

We search then for such a modification of (5) as will provide stability with mini- 
mum degradation of accuracy. The following discussion will establish that a usable 
and in fact essentially optimum modification of (5) consists of replacing the series 
of six coefficients 1/6, 1, 15/6, 20/6, 15/6, 1 multiplying the [ ] by new constant 
coefficients Y = 95/288,1, A = 25/24, B = 35/72,C = 5/48, D = 1/120 respectively 
and leaving (5) otherwise unaltered. It is interesting to note that the ratios of the 
new coefficients to the old form a rather strongly decreasing sequence: 1.98, 1, 0.42, 
0.15, 0.042, 0.0083, which reminds one of the well known technique for stabilizing 
electrical filters involving feedback by somewhat enhancing the low frequency gain 
and strongly depressing the high frequency gain. 

In searching for an appropriate modification of (5) it is inadvisable to tamper 
with the coefficients not pertaining to the [ |], and this will be borne out by later 
analysis, for these coefficients are clearly just such as to make the integration of a 
5th degree polynomial y(x) come out exact (the [ | will vanish for y(z) a 5th degree 
polynomial). However, the coefficients multiplying the [ ] have no such unique 
significance and we are free to modify them to suit our purpose. 

To dispose of the possibility of generalizing the coefficient 1 in the second equa- 
tion of (5): So long as this coefficient remains 1 we can delete the second equation 
entirely from the considerations as being merely an identity, and we ultimately do 
just that. In the interests of generality the author has experimented some with 
modifying this particular coefficient numerically and has indeed found that any 
value other than 1 for it, beside costing an additional multiplication, degrades both 
the accuracy and the stability. 

The remaining 5 equations of (5) with the coefficients 1/6, 15/6, --- 1 replaced 
by arbitrary constants Y, A, B, C, D, may then be studied for stability by introduc- 
ing a small variation of each of the 5 independent quantities (y, ha, hb, he, hd), 
namely (dy, dha, dhb, dhe, hd), and studying how this latter quintuple changes as 
we integrate from x to x + h [3]. The quantity f is to be regarded as not independent 
but a function of y in virtue of the differential equation. After some calculation we 
find that the quintuple (éy, héa, héb --- hid), regarded as a 5-component vector 
V(x), obeys the equation 


(6) Viz +h) = MV(z) 
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where M isa 5 X 5 matrix: 



































M = 
P ¥(p — 2) Y¥(p — 3) ¥(p — 4) Y(p — 5)| 
1 1 4 1 gt — 
ees; ot “Sie wa ¥> oe ie 
_Ap A@-2 ,, A@—3) Alp ~ 4) A(p ~ 5) 
1—Yp AT 2 ir PT ot ore eh es 
Bp* B(p — 2) B(p — 3) B(p — 4) B(p — 5) 
@) 1—Yp 1—Yp at E> nie OE eT 
Cr C(p — 2) C(p — 3) 1402-9 5, ce-5) 
1-—Yp 1~—Yp 1— Yp 1—Yp 1—Yp 
Dp? D(p — 2) D(p — 3) D(p — 4) D(p — 5)} 
1-—Yp 1-Yp 1—Yp 1—Yp 1— Yp | 
with 
af(x, y) 
= pry? 
(8) p ay 


We note that the 5-dimensional vector space of V and M is a different space from 
the n-dimensional space of y, f, a, etc. 

We have treated p as though it were a scalar quantity even though for n > 1 
it is really ann X n matrix h(df;/dy;); but it is only the smallness of p, insurable by 
appropriate choice of h, which is important in our argument, not its matrix charac- 
ter. The difference between p(x + A) and p(x) has also been neglected, for it gives 
rise to errors involving one factor h more than we need consider. 

The characteristic equation 0 = | \é,, — M,, | of M turns out to be: 


0 = (1 — Yp)(A — 1)° 
+ (24+ 3B 4+4C0+5D-—-(14+A+B+C + D)pj\(a — 1)‘ 
(9) + (6B + 24C + 70D — (2A + 6B + 14C + 30D)p}(A — 1)° 
+ [24C + 180D — (6B + 36C + 150D)p](A — 1)° 
+ [120D — (24C + 240D)p](A — 1) — (120D)>p. 


One root of this equation, which may be found by substituting a power series in p 
into it, and which we shall call the principal root Xo , is essentially a function of p 
only, depending but slightly on Y, A, B, C and D: 


6Y —3+A — (1/5)Cp’ 
D 6! 
—49 + 105Y + 144A 7/5)B — (14/5)C — Dp’ 
+ 9) + si /5) a how F + 0(r'). 





No =e? + 
(10) 
+ 





This is a consequence of retaining the coefficients in (5) not pertaining to the [_ |. 
The root Xo is thus essentially a property of the differential equation system (1), 
and whether or not it lies inside the unit circle in the complex \ plane determines 
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whether the subject system, as distinguished from our numerical method, is stable 
or not. On the other hand, the four further roots of (9), which we shall call “ex- 
traneous” roots, depend strongly on A, B, C, D and only weakly on p and Y; their 
location relative to the unit circle determines the stability of the integration method 
itself. These roots must lie inside the unit circle for stability of the method, and the 
nearer they are to the origin the more stable the method will be. 


4. Determination of Parameters. The parameters Y --- D are now to be chosen, 
primarily to optimize the stability of the method and secondarily, if any freedom is 
left over, to optimize the accuracy within the restriction of optimum stability. The 
author regards optimum stability as essential to an automatic general-purpose 
method, for the rapid elimination of disturbances characteristic of good stability not 
only makes an automatic starting process feasible and permits accurate integration 
across finite discontinuities of f, as we shall see below, but also minimizes the error 
due to interaction of disturbances with non-linearities of the differential equations.t 
Since there are four extraneous eigenvalues whose locations in the complex plane we 
wish to control and we have five parameters free, we can expect to have considerable 
control over stability and accuracy. What actually happens is that A, B, C, D 
determine stability and Y is left free to optimize accuracy. Thus we can arrange for 
a truncation error of 0(h') even though we are using 5th degree polynomials, the 
explanation being that in each integration step we use both the 5th degree poly- 
nomial available at the beginning and the one available at the end of the step. 

Now it is easy to bound | p| (bound the magnitudes of its eigenvalues if it isa 
matrix) by control of h during the numerical integration process, while it is much 
more difficult actually to compute p for n > 1. Therefore it seems best and is cer- 
tainly simplest to choose Y, A, B, C, D independent of p, i.e. as absolute constants, 
in such a way that stability is guaranteed for as large a range of p as possible. This 
is substantially accomplished by considering (9) with p = 0 (whereupon Y drops 
out, indicating that it has little influence on the stability of the method) and then 
choosing A, B, C, D so that the four extraneous rocts coincide at 0. Thus, we require 
(9) for p = 0 to take the form (A — 1)\* = 0, and it does that for A = 25/24, 
B = 35/72, C = 5/48, D = 1/120. The choice of Y is then made to nullify the 
coefficient of p*® in (10), which has no effect on the stability but optimizes the 
accuracy. This determines Y = 95/288. For stability for p = 0 we then depend on 
the fact that the extraneous roots are continuous functions of p, so that they cannot 
move very far from the origin provided p is appropriately limited. 

In order to get a better picture of the behavior of the extraneous roots as func- 
tions of p, we first note that for small p they are the roots of 

: 3 
(11) is Teo ? 


as can be read off from (9) with the chosen values of the parameters inserted. It is 


+t Areport by E. Fehlberg [4], has just come to the author’s attention. Fehlberg exhibits 
other choices of parameters which produce smaller truncation error than Adams’ and the 
author’s choice, but at the expense of much poorer stability, cf. Fehlberg’s tables 3 and 4. For 
m = 5 the gain in computing speed for the same error is greatest and is (1/0.0801)'” = 1.43, 
which the author considers not worth the risks incurred with the much poorer stability. 








30 ARNOLD NORDSIECK 


fortunate that the numerical coefficient in (11) is so small, for the p'* deycz ence 
of the roots is a rather strong dependence. The roots have also been computed for 
p a real number between —1 and +1, and these are shown in Figure 1. We see that 
stability will be guaranteed with a comfortable margin of safety if the interval is so 
chosen that p lies effectively inside the dashed curve. This boundary corresponds 
to | Yp| < 1/8, which is a convenient form of test for a computer. 

The author has done considerable searching for other favorable choices of A, B, 
C, D with the thought in mind that if the extraneous roots never coincided they 
might move away from the origin more slowly as | p| increased, than they do 
according to (11). However, all other choices tried were inferior in point of both 
stability and accuracy. 

The choice of parameters made above seems optimum among choices restricted 
to constants independent of p. The potential advantage of a more elaborate pro- 
cedure in which the matrix p is numerically computed at every stepand Y,A --- ,D 
are made chosen functions of p, implying a nonlinear process tailored to the subject 
differential equation system, is an interesting topic for future investigation, for it 
might lead to faster (though less accurate) methods of solving some classes of 
equations. 

The working equations of the method have now been determined completely and 
they are summarized in Appendix A, equations (A4). 

The working equations having been determined, the precise connection with 


|Al=! 








\ 


Fig. 1.—The extraneous roots of the characteristic equation as functions of p for real p, 
plotted in the complex \ plane. As p departs from zero these roots depart from the origin along 
the loci shown. Loci marked + correspond to positive p and loci marked — to negative p. 
Counting outward from the origin along each locus, the points plotted represent in order, 
|p| = vs, 4, 4, 3, 1. The real positive extraneous root coalesces with the principal root at p = 
—0.88, producing a conjugate pair. The dashed curve encloses all extraneous root values per- 
mitted by the interval control tests, which limit |p| to values < .38. 
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other methods can be deduced by ascertaining the equivalent quadrature formula 
for the method. This can be done by expressing the part h{ } of the first working 
equation in terms of past values f(z — h), f(z — 2h), ete., by repeated application 
of the working equations. We find that the equivalent quadrature formula is 


h 
(12) y(z +h) — y(z) = 1440 {475f(z + h) + 1427f(x) — 798f(x — h) 


+ 482f(2 — 2h) — 173f(x — 3h) + 27f(x — 4h)} 


which agrees exactly with the Adams formula of corresponding degree. By way of 
confirmation of this conclusion we observe that the characteristic equation for 
small variations in the Adams method coincides with (9) when the chosen values 
of the parameters are inserted into (9). 


5. Change of Interval. We indicate how to perform the three useful changes of 
interval: h’ = —h,h’ = Bhand h’ = 6h (where in binary computer operations 8 
is preferably taken equal to 2): 


Reversal Increase Decrease 

—h Bh Boh replaces h 

y y y PRO 

f f f pie 

(13) —a Ba Ba * a 
b B*b Bb “ b 

—c Bec B-~*c ™ c 

d Bd Bd “ d 


The rules for changing a, b, c, d are clear from (3). 

The simplicity of the rules for changing the interval is evident here. 

Every change of interval of any of the three types induces a disturbance in the 
system, but the disturbance affects mainly the higher derivatives and clears out in 
a few steps because of the choice of parameters. These transient phenomena will be 
described in more detail in the next following section. 


6. Behavior of a, b, c, d. A qualitative understanding of the behavior of the 
quantities constituting the method’s “memory” is required in order correctly to 
design the interval control logic and the starting procedure. 

We first describe the ‘‘normal’”’ or steady behavior which prevails when no 
transients have been induced by interval change or f-discontinuity or otherwise, 
within the preceding 4 to 8 steps or so. Then the quantities a, b, c, d ‘“‘lag”’ behind 
the current value of x, a a little, b more, c still more, and d most, in the sense that 
they equal the ‘‘true” higher derivatives of y evaluated at points z — 6h, where 
0 < @ XS 2. This lagging behavior is related to, and is in fact a necessary consequence 
of stability. A close analogy exists between this and the “‘stable physically realizable 
filter’? of electrical engineering theory, and likewise the causality discussions in 
physics. The indicated behavior may be established (and incidentally some formulas 
of later use for deriving the truncation error found) by assuming that a 7th degree 
polynomial y = P;(2x) satisfies the differential equation exactly and that f, a, 6, c, d 
are corresponding polynomials of 6th --- 2nd degree, and solving the working 
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equations (A4) for the coefficients by some rather lengthy algebra. The result is 
(2) = Aya) — 72% ya) + 840% (2) 

a(z) = oT y" (xz) — 61 y (zt) + 7 y (2 

hr? mt hr 
31 y” (x) — 100 éi y 
h° sere h° 
c(x) = nY (x) — 52 1/2 a 


6 
b(x) = "(z¢) + 1110 2/3 i y (x) 
(14) ¥ 
VI ror h Vil ‘ 
y (x) + 525 = y""(2) 


h* wy h° 
d(x) = 5 y’""(x) — 12 rT y (x) + 91 7 y (2). 


These formulas are then in error by 0(h’) for any general y which is differenti- 
able sufficiently many times. The last of the four formulas shows that d(x) = 
i y’”"" (x — 2h) + 0(h*), so that d lags by very nearly two steps. 

We may describe this ‘normal’ behavior in another way, namely, by observing 
that the polynomial evaluated at z is always essentially the polynomial fitted to 
the values of y at x, x — h,x — 2h, x — 3h,x — 4h. The 5th derivative of this 
polynomial naturally agrees best with the 5th derivative of the true solution y at the 
mid-point of the fitting interval, which explains the last equation of (14). The close 
relation of our method to the Adams method also becomes clear from this point of 
view. When we advance from x to x + h the working equations in effect change the 
old polynomial fitted at x — 4h --- x into one fitted at x — 3h --- x + h. In the 
approximation that the 4th powers of the extraneous characteristic roots may be 
neglected, all disturbances clear out in precisely 4 steps, corresponding to the 
memory of the method having a “‘time-span” of just 4 steps. Thus, we have ar- 
ranged effectively to keep and use what Adams actually keeps and uses, namely the 
last four previous ordinates, whereas actually we keep quantities much more suit- 
able for interval modification. 

As for “abnormal” behavior of the remembered quantities, the simplest im- 
portant case of this occurs upon reversal. The quantities exhibit a hysteresis after 
reversal, most pronounced in the case of d(x) which has the most lag. The behavior 
of d in reversal is illustrated in Figure 2, which shows essentially that d stays quite 
strictly constant for four steps after a reversal, then abruptly resumes normal be- 
havior. Since what was a backward-fitted polynomial before reversal becomes a 
forward-fitted polynomial after reversal, we may say that d and indeed the poly- 
nomial as a whole “‘freezes’”’, remains the same and marks time until enough steps 
have been executed for it to become normal for the current point x, then behaves 
normally. 

The other important type of abnormal behavior is the response to shock excita- 
tion. Shock excitation occurs severely in starting, when the normal a - -- d are not 
known; mildly enough to be harmless in increasing or decreasing the interval, when 
the main terms but not the “lag” terms in (14) are correctly modified by the simple 
rules (13); and more or less severely when f has a discontinuity so that the change in 
the polynomial is large in one step. Here again d(x) shows the most violent behavior 
and its behavior in all shock-excited transients is essentially an oscillation lasting 
just four steps. 
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TABLE 1 

: | y | f | 24a | 72b | 48c 120d 

te 0 ERS Se os 0 
ee a 2a ras hd | A | 25A | 35a | 5A A 
xo + 2h | 3+ rae ha | 4 | -23a | ~69a | ~130 | —3a 
to + 3h 2 vata Ja | 4 | +138 | +450) na) 3 
tet th | ¢ + i404 | 4 | -30 | -a | -3a| a 
zo + 5h | Sha SY te he 0 | 0 





In order to become familiar with the detailed behavior of such a transient, we 
treat a simple case which approximates the general case of an isolated discontinuity 
of f with finite jump: Let y = f = 0 for z S 2 and assume that a --- d have their 
normal values of zero for x S x .Letf = A= constant for z > 2 and apply the 
working equations (14) five times in succession and Table 1 results. 

! ! 
Evidently = a, s b, e c and iid are behaving like numerical approximations to the 
“§-function” of x and its first, second and third derivatives respectively. Meanwhile 
the transient in y, represented by the terms with denominator 1440, is a decreasing 
oscillation also lasting just four steps, and the ultimate value of y is exactly what 
one would get by connecting the last point sampled at which f = 0 with the first 
point sampled at which f = A by a straight line segment. This essentially best per- 
formance in integrating across a discontinuity is unique to our choice of parameters. 
A reasonable upper bound for the magnitude of the error in y due to such a dis- 
continuity is 4 | hA| where A is the jump in f. One can hardly do better without 
sampling in between these two points, i.e. decreasing h; but by controlling h one can 
bound this error. 


7. Estimation of Errors. In discussing errors in the solution y(x) we must distin- 
guish between the error present at the beginning of an elementary step and the 
error contributed by the execution of that step. The error present at the beginning 
of a step, sometimes called inherited error, is the net result of all the errors con- 
tributed by all the previous steps, each modified according to the action of the 
differential equation between the point of origin 2x’ of the error and the current 
point x. Letting E(x’) represent the error contributed by a step of length h taken 
at x’, we can write for the inherited error at z if the integration began at 2 : 


(15) E(x) = 0 E(x’) TI 02”) 

z'=z9 2’ '=z' 
where \» & e” is the principal root (10). The sum involves the summand once for 
every elementary step taken and similarly the product. Equation (15) illuminates 
the relation between inherited error and error contributed by an individual step. 
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The product in (15) may also be written in adequate approximation as: 


(16) II do(z”) = exp | of (x”) ax"\ 
a 2’ \ z’ oy 

which shows that the product is a property of the differential equation, independent 
of integration method and of interval choice. It is clear then that although by careful 
design of the method and choice of interval we may be able to reduce E(x’) down 
to about half the least count in the register (but no further because of inevitable 
rounding), nevertheless such measures have no effect on (16). Consequently if A 
is the largest eigenvalue of the matrix (16) the error at the conclusion of the inte- 
gration will be in general at least about 4 | A | times the least count. The number of 
correct significant digits may at most be preserved through the calculation if the 
magnitude of the solution increases by A or more; if not, the significance (i.e. the 
number of correct significant digits) will decrease. If a problem has A > 8”, where 
L is the number of base 6 digits in the register, then it is useless to attempt the 
problem at all by fixed point arithmetic, for there will be no correct significant 
digits left at the end of the calculation. Floating point could help if the magnitude 
of the solution increases meantime; if not, nothing wil! help except increased 
register length. 

We have dwelt on the above points because they show that the best that can 
be done with any method is approximately to preserve the number of correct sig- 
nificant digits in the solution, and this essentially defines a best or optimum method. 
Some of the test examples exhibited in Section 10 below show nearly complete 
preservation of significance through as many as 10° steps and with A as large as 
10° or so. 

Turning now to discussion of E(2), we assert that the contributions to H(z) are: 
a) truncation error incurred by terminating the formulas (A1) to (A5) with a given 
power of h; b) discontinuity error incurred in integrating past a discontinuity of f 
(ef. Section 6); c) iteration error resulting from incomplete iterative solution of 
the implicit equations for y(z + h); and d) roundoff error resulting from using 
registers of finite length to perform the arithmetic. 

The truncation error may be found by making the same assumptions y = P;(z) 
etc. as were made in deriving equations (14) and calculating y(z + h) — P7(x + h) 
— y(x) + P;(x), using the first and second working equations and (14). We find 
that the truncation part of FE, which we call F, , is given by: 


h’ 
(17) E,(a + h/2) = 72 = y (x) + O(h’). 

It is interesting to note that the truncation error is closely related to the principal 
root of the stability matrix. In fact, if we replace p arbitrarily by the operator 
bd (because the proof that p is precisely equivalent to h $ is not apparent), then 
e” becomes the “true” displacement operator ¢ and o(p) becomes the approxi- 
mate displacement operator of the method. Thus the difference \o(p) — e” with p re- 


h(d/dz) 


placed by ho , and applied to y(x), would seem to yield the truncation error. The 


term in f° in the truncation error was determined by exploiting this relationship, 
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yielding: 

/é "76 h’ vil h° vill 9 
(18) E,(z + h/2) = 25 y (zs) - 440 5 y (x) + 0(h) 
The discontinuity error, called Ez , is bounded by the inequality 
(19) | Za(x)| S 3 | A(f(r+) — f(z_))| 


as we saw in Section 6. 

The iteration error, called £; , depends on how we solve the implicit equation 
system, and we choose to solve it by doing just two iterations, or more precisely: 
We calculate a first trial value y (xz + h) from equation (1) of (A4) with the 
[ ] term left off (the “predicted” y(z + Ah) in Milne’s terminology); calculate 
f(a + h) = f(x + h, y (x + h)) and insert it on the right of the complete 
equation (1) of (A4) to give an improved y®(z + h); and repeat the procedure 
just once more, so that by definition in this method the final values of y(z + h) 
and f(z + h) are y® (x + h), respectively f(x + h, y®(x + h)). The reasons for 
choosing so are that f need be calculated only twice, that the convergence of the 
iterative procedure and the (related) bounds on p can be estimated from two 
iterations but not from less than two, and that the iteration error is sufficiently 
small. For the special case n = 1 (a single first-order differential equation) one can 
do better by solving the implicit system by interpolative methods with the same 
number of computations of the derivative; for general n, however, one would have 
to compute the derivatives 2n times at least in order to apply interpolative methods, 
which we regard as uneconomical. The convergence is determined by the equation 
re (3) (2 (2) a » 
(20) gf mig a Fpl? — 9); Fie 


and the “iteration error” in y(z + h) by 

(21) E; = y® — y® = (¥p)(y™ — y°) & —Y'p*h'y(z) 

which is proportional to h* with a small coefficient so long as | Yp| < } as we shall 
require, and is therefore overshadowed in general by the truncation error. Equations 
(20) and (21) follow from iterative treatment of equation 1 of (A4). 

The roundoff error E, , finally, is determined by the care with which both the 
computation of derivatives and the computations of (A4) are done, and with 
sufficient care can be as small asabout } the least count in the effective register used 
and approximately statistically independent from step to step. The author has 
found it best to keep logs(| h |’) extra “guard” digits in y, above and beyond the 
number kept in f, a,--- d, in order to minimize the accumulation of roundoff 
errors in y when the number of elementary steps is large. 


8. Automatic Interval Control Logic. In order to describe the interval control we 
must first outline the 3 stages in which a step x — x + h is performed. Stage 1 con- 
sists of “predicting” all six quantities y, f, --- d at zx + h, ie. applying equations 
(A4) without the [ ] terms, using a tentative value of h. The first tentative value 
of h actually tried is the value which was accepted in the last previous step or the 
next larger value if the conditions (given below) for increasing h were fulfilled. Note 
that Stage 1 is exactly reversible in a digital machine, so that if h later turns out to 
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be wrong the beginning values of y --- d can be exactly recovered without the 
need for additional registers for saving them. Stage 2 consists of solving the implicit 
equation system for y(z + h) and f(x + h) by iterating twice as explained in the 
preceding section. This stage is not exactly reversible and 2n registers are therefore 
provided for saving the beginning values of y and f. At the conclusion of Stage 2 
enough information has been developed to decide whether the interval tentatively 
being used is small enough; if it turns out to be not small enough the beginning 
values of y --- dare recovered, the interval is reduced (by a factor 8’ = } ina 
binary computer) and Stage 1 is again entered. If the tentative interval is found 
adequate we proceed to Stage 3, which consists of “correcting” a, b, c, d by adding 
the [ ] terms. 

Two testsare made at the conclusion of Stage 2 and failure of either signifies that 
h is too large; the two tests are respectively 


(22a) | yi” toe ye? _= ~* $ | ys” TTL ys” | max 
and 
(22b) | fix + h) cs f? — s B*/| h | 


where e is a specifiable positive integer and ‘‘max’’ means the largest of the n 
components i = 1, 2, --- n. It is clear that these tests are first possible at the 
end of Stage 2, since they involve quantities developed only in that stage. While the 
tests are being made it is also determined whether both tests are “‘over-satisfied”’, 
i.e. so well satisfied that the next larger h would likely also satisfy them, and if so 
the interval may be tentatively increased for the next following step. 

Satisfying test (22a) insures that the largest eigenvalue of p does not exceed 
0.38 in magnitude (cf. equation (20)) and, therefore, that the stability is good 
(ef. Figure 1) and also that the iteration error is small enough to be overshadowed 
by the truncation error (cf. equation (21) ). The test is not formulated in the ideal 
way, which would be to require the Euclidean norm of the difference vector to 
decrease by 4; instead we require that the largest component of the difference vector 
decrease by at least $, which is equally effective in insuring convergence, works for 
any order n, and requires less computation and less registers. 

Satisfying test (22b) then has the effect of roughly bounding the truncation 
error and the discontinuity error in such a way that the accumulated error in inte- 
grating a standard distance Ar (which we take equal to 1) is independent of the 
elementary step-lengths used and about equal to 6 “. In effect, instead of having to 
specify the elementary step-lengths to be used, the programmer tells the com- 
puter he wants the eth digit in y to be correct after integrating a unit distance along 
the z axis and the computer is expected to choose the elementary intervals to achieve 
this result most economically. Note, however, that in this connection the discussion 
of preservation of significance for unstable equations given at the beginning of 
Section 7 must be kept in mind. 

Test (22b) is derived from equation (17) by the following rough argument. 
We divide the interval (xo , 7» + 1) into subintervals in such a way that within each 
subinterval h is constant. Then summing (17) over the kth subinterval gives: 
h® Zk+1 

& 


(23) & = 3 Bj (2) ee 
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Now the computation provides an estimate of h°y”’, namely, h[f(z +h) — f’}, 
as may be deduced from the 6th equation of (A4). We use this estimate to bound 
hy’ for all x by requiring satisfaction of test (22b) in every elementary step. 
Thus h* | y"’ | < 6‘ and the accumulated error is, roughly speaking, bounded by 





24 no. of subintervals 
(24) | » &| <x 0 


which is not likely to be much greater than 6 °. We see also that the general effect 
of bounding h°y”’ is to cause each part of the total integration interval to con- 
tribute to the error in proportion to its length, which tends to minimize the total 
number of steps to achieve a given accumulated error. The argument is necessarily 
somewhat crude, for we cannot do what one would ideally like to do, namely, 
bound h*y’”’, because there is no estimate of it available (without increasing the 
degree of the method). Test (22b) also bounds the discontinuity error, equation 
(19), for a discontinuity if f clearly appears directly in [f — f”], so that bounding 
h{f — f”| just bounds (19). 

In addition to availability of an estimate there is a further practical reason for 
formulating test (22b) in just the way shown, at least in a fixed point arithmetic 
operation, namely, that it permits the widest possible range of choices of h without 
either member of the inequality falling outside register range. If one wants to inte- 
grate across large discontinuities of f and still be free to demand accuracy of the 
order of the least count, it is clear from (19) that h must be reducible to or near the 
least count; on the other hand, for maximum size steps when f varies slowly and 
smoothly Ah must be increasable to or near the greatest count of the register. In 
practice the author has had the interval vary all the way from 2° to 2™ in a 39 
binary digit machine. 

In the main then, the interval is selected by requiring it to be the largest interval 
satisfying both tests (22a) and (22b). However, four minor modifications of this 
basic rule are introduced in order to improve the usefulness and efficiency of the 
method and the smoothness of the automatic interval control, as follows: 

Since the programmer cannot predict what intervals will be used he is given the 
privilege of specifying a maximum interval hy , so that he has assurance that the 
solution will be available at least at the points zx» + (integer)ho . The automatic 
interval control then includes a feature preventing an increase in the interval 
whenever such increase would result in skipping over one of the above points zx» + 
(integer ) ho . 

Next, when any considerable amplitude of shock excitation has occurred it seems 
best, judging from Table 1, to choose the interval at the onset of the shock, then 
leave it unchanged until the transient due to the shock has subsided. In fact, if the 
interval is changed while the strong transient is still present this interval change 
itself results in a new shock excitation, and the interval control tends to become 
erratic in the sense that the interval is reduced too much and for too long, a phe- 
nomenon which the author has observed experimentally. The interval control itself 
contains feedback loops, we may say, which can cause erratic behavior, although 
not genuine instability because the computer takes refuge in reducing the interval 
in response to any uncomfortably large disturbance. The main rule, if not modified, 
leads to just such behavior because, as we see from the last column of Table 1 


-p° 
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the change in d is 4 to 6 times greater in steps subsequent to the first step after 
onset of the disturbance than in the first. To avoid this misbehavior the computer 
is programmed to recognize the characteristic A, —4A, +6A, —4A --- pattern 
and to leave the interval unchanged on the 2nd, 3rd and 4th steps provided they 
conform to this pattern within certain tolerances. This effectively prevents the 
interval control from interfering with the expeditious elimination of transients and 
results in preserving the ideal accuracy and speed represented by Table 1. 

Another form of undesirable interference from interval control occurs in connec- 
tion with reversal. Suppose that reversal has just occurred and that test (22b) is 
dominant in determining the interval, as it often will be. From Figure 2 we see that 
just after reversal d stays constant for 4 steps. This means that (22b) will be over- 
satisfied and the interval will be increased, whereas it should clearly not be increased 
since we are retracing steps for which the interval was presumably already correctly 
chosen earlier. The subsequent behavior would involve an unusually large shock 
when the “slack’”’ in d is eventually “taken up” and an unnecessarily large interval 
decrease, again a phenomenon the author has observed in practice. The remedy for 
this misbehavior is simple: we program in a rule preventing interval increase for 
the first four steps after any reversal. 

Finally a rather interesting type of misbehavior can occur when f tends toward a 
constant or indeed toward any 4th degree polynomial after an earlier more violent 
behavior which required a small interval. In these circumstances we want and 
expect the interval to increase rapidly, but if the parameter 6 “ of test (22b) is very 
small, say only a few times the least count, then such increase may be prevented 
entirely by persistent roundoff noise in the “remembered” quantities. If f tends 
asymptotically to a 4th degree polynomial d should tend to a constant and 
[f(z + h) — f”| should tend to 0. What happens then is that so long as roundoff noise 
persists either (22b) is barely satisfied and the interval is not increased, or if (22b) 
is oversatisfied and an interval increase is attempted the roundoff noise in d is 
magnified by a factor 6* according to (13) and causes test (22b) to fail on the next 
step. Now, unless special measures are taken, the roundoff noise can indeed persist 
and prevent interval increase indefinitely. Thus we may get into (and the author 
has actually got into) the absurd situation of taking 4000 steps to integrate * =0 
from x = } tox = 1 (provided f was non-zero for x < 4). The remedy for this mis- 
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Fig. 2.—Hysteresis behavior of the “‘remembered”’ quantity d. The dashed curve is the true 
value ofd(z) ; the solid curves show the behavior of d in the computation. 
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behavior is not modification of the rules for interval choice, but a peculiar, carefully 
chosen rounding procedure for the multiplications by A, --- D involved in the 
working equations so as to guarantee that roundoff noise will disappear in a finite 
(and minimum ) number of steps just as other transients must and do because of the 
stability. The discussion of choice of rounding is rather long and is also of interest 
for other iterative procedures in numerical analysis, therefore, it is given sepa- 
rately in Appendix C. 

The main rule amended by the four modifications just described provides a 
stable, non-erratic and generally reasonable behavior of the interval size in all 
cases which have been investigated, and the cases investigated were purposely 
chosen extremes in which the interval had to vary rapidly and widely. The interval 
still does not increase as fast when it should increase as it decreases when it should 
decrease, but this is hardly avoidable since both the finite rate of clearing of tran- 
sients and the requirement of not skipping over the points z» + (integer)he act to 
delay interval increase. 

If, when h has been reduced to the least count, test (22b) still fails, a programmed 
stop is encountered. Almost any major malfunction of program such as overflow 
in the computation of derivatives or elsewhere leads quite immediately to this 
programmed stop because of the extreme sensitivity of test (22b). 


9. Automatic Starting. The essential idea which makesautomatic starting feasible 
is that if we set off with entirely abnormal values of a, b, c, d, say putting 0 for each 
of them in the absence of any evidence as to their normal initial values, then upon 
integrating a few steps they will assume approximately their normal values if the 
stability is sufficiently good. Such a method of starting has the advantage of using 
mostly the normal integrating program, which has to be supplied in any case, 
requires very little extra programming of special nature, is of use only during 
starting. Since a modern computer can execute at least about one step per second 
in even rather complicated differential equation problems, the start can be accom- 
plished blunder-free and accurately in a matter of seconds or at most minutes. 

Several complications must be dealt with in providing a satisfactory automatic 
start: the proper interval for the first step forward from x» is not known in advance 
any more than are a, b, c, d. There is a certain degree of incompatibility between 
automatic starting and interval changing since the starting essentially involves 
eliminating a very large transient and, as we saw in the preceding section, changing 
the interval during a large transient can lead to erratic interval behavior. In any 
case, application of test (22b) during the first few steps of the starting process 
would be meaningless since the test was derived on the assumption that a --- d 
had nearly normal values; this is illustrated by the fact that when a, b, c, d are zero 
the quantity [f(z + h) — f”] is O(A), not 0(h°) as it normally is. Finally, although 
one would ideally like to use points to the left of zo for starting, corresponding to 
fitting a polynomial to the left of z> and thus obtaining what we have called normally 
lagging values of a(zo) --- d(xo), this cannot be done because it would imply that 
f is defined to the left of xo , as it may not be. 

The detailed schedule of the starting procedure will now be described and in the 
process the way in which the complications listed above are dealt with will become 
clear. The overall objective of the starting procedure is to fit a 5th degree poly- 
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nomial for y to the points ro , %» + h/2, ro + h, xo + 3h/2, xo + 2h, thus determin- 


ing a(2xo), b(ao) --- d(ao), where h is the correct interval (also to be determined ) 
for the first step xo — 2 + h. 
First we set the initial values y(zo) = y’ aside for safekeeping, set a --- d 


equal to zero and do a tentative step forward x» — 2» + ho , where ho is the maximum 
interval permitted. Test (22a) (but not (22b)) may now be applied since its opera- 
tion is essentially independent of whether a --- d have their normal values. If 
(22a) fails the interval is reduced, the beginning values at 2» are recovered and a 
shorter tentative step forward from 2 is taken, the program used here being just 
the same as in normal integration. This process continues until an h has been found 
which satisfies (22a). 

When (22a) has been satisfied three more steps forward are taken, followed by 
a reversal and four steps back to 2» , all eight steps being taken at a constant inter- 
val. The reason for taking just four steps either way is that it provides just enough 
information to determine a 5th-degree polynomial. 

We are now back at x» with a value of y somewhat in error but with first approxi- 
mations for a --- d which are already good to a fraction of a percent because of the 
high degree of stability of the method. The correct value of y(zo) is reinserted, the 
sense of integration again changed to forward and another four steps forward and 
four steps back to z» are taken, all at the same constant interval. 

During the last backward step listed (the 16th step of the starting process) test 
(22b) is activated, for now the quantities a --- d are so nearly normal that this 
test is significant. Test (22b) must be made neither too early during the starting 
process, for then [f(z + h) — f”] is not yet 0(h°); nor too late, for as the process of 
integrating four steps back and forth is continued, [f(z + h) — f”] tends to zero 
in any case (refer to the hysteresis behavior of d described in Section 6). Thus 
there is a sort of psychological moment for doing test (22b) during the starting 
process. The author has found by “experimental mathematics” that [f — f”] is 2 to3 
times larger on the 16th starting step, for all equations and all h’s, than it is in the 
ultimate normal integration process. Thus, applying the test at this point results in 
a slightly conservative initial choice of h. 

If (22b) is not satisfied the interval is reduced and we go back to the very begin- 
ning of the starting process. If (22b) is satisfied, y° is reinserted, the sense of inte- 
gration is changed to forward and the starting process may be considered almost 
completed. In fact, for all cases except those with very high accuracy requirements 
and very unstable equations the above process provides a satisfactory start. In the 
exceptional cases one can do a little better, typically a factor of six in the initial 
truncation error, by extending the starting schedule to include four more steps 
forward and four back at half the interval eventually to be used (making 24 starting 
steps in all after both tests are satisfied) and we actually take these extra eight steps 
in order to be quite sure that errors attributable to starting are less than the normal 
running truncation error. More precisely, after test (22b) is satisfied during the 
16th starting step, we reinsert y°, change the sense to forward, halve h, integrate 
forward four steps, reverse, integrate back four steps to zo , reinsert y°, change the 
sense to forward, double h and now regard the starting process as complete. 

The chief effect of performing the last eight starting steps at a reduced interval 
is to reduce the amount of lead in a, b, c, d, which is beneficial because they should 
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lag, as they would if we had used a backward fitted polynomial. In any event, the 
truncation error in the first step after starting as above is less than normal. 

Note that we have avoided ill effects due to changing interval during a transient 
by insisting that if any starting step at all is taken with a given interval h, at least 
eight are taken without changing A. 

The transients during the early stages of the starting process are often large 
enough to cause overflow of the computer registers, and it is interesting to observe 
that such overflow will do no harm, for test (22b) is a very sensitive test and will 
almost certainly be violated if there are any previously occurring overflow errors. 
When this test is violated the computer simply discards all its previous computa- 
tions, including any overflow errors, and starts afresh with reduced interval. The 
author has observed this effect many times, always without ultimate consequences. 
Persistent overflow caused by incorrect scaling of x, y or f is of course another mat- 
ter, but one which comes to light very quickly in the form of the programmed stop 
mentioned earlier. 


10. Test Problems Done by This Method. The differential equation problems 
used to develop the program and to rectify programming errors were those for the 
sine function and the exponential function. The normal truncation error for these 
‘“‘well-behaved” problems was found to agree with (18). 

A test problem to exercise the automatic variable interval feature thoroughly 
and to verify the behavior for discontinuous f was then devised as follows: 


4 (0 for |x —}3|22™ 


os) dx \2° for |x-—}3|<2™ 


to be integrated from 0 to 1 with ho specified as 2° and 8“ specified as 2". This 
involves having the computer search the z-axis efficiently for an extremely narrow 
region in which f ~ 0, finding the area under the curve in this narrow region very 

















TABLE 2 
Steps | h x | y-2* | “Correct” y-2” 
0 } 2s | 0 | 0 | 0 
157 as | 1/2-2" | 0 0 
169 | 2s 1/2 015 564 | = .015 564 
176 | 2 | 1/242" | .031 149 | .031 128 
370 by: Sve of 1 | .031 128 | .031 128 
TABLE 3 
Steps h x | ya Error 
| | — 
0 2-8 | —1/2 | 0 | 
202 2-3 | —2-s .098 177 .000 002 
214 | 0 | 196 352 .000 003 
227 a2 2% | 294527 | .000 003 
505 2-8 1/2 .392 700 .000 001 
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accurately and then searching the rest of the z-interval at high speed again. The 
performance on problem (25) is shown in Table 2. 
The “‘correct” y means the exact area of the figure obtained by joining the consecu- 
tive pair of points sampled, with h = 2”, at which f changes, by a straight line. 
The interval actually increased 64-fold temporarily between corner and center of 
the curve. The somewhat slower recovery of the interval on the increasing-interval 
side is exhibited in the difference between 194 steps from x = 1/2 +2 tox = 1 
and 157 steps from x = 0 tox = 1/2 — 2. The recovery of the interval to hy = 2° 
at all is evidence that roundoff noise does not persist in the ‘cremembered’’ quantities. 
A test problem similar to the above with a very narrow but smocth analytic 
curve was also treated: 


(26) yg 2 : 

dx x? + (2-”)? 
to be integrated from x = —1/2 tox = +1/2 with ho specified as 2° and B~* 
specified as 2~”. The result of this computation is given in Table 3. 

The interval evidently did not have to decrease so much in this case because of 
the smoother curve to be integrated. The same comments in regard to increasing h 
apply here as in the previous example. The accumulated error is much less than 
2” because of the simple symmetrical character of the curve being integrated. 

Next, a typical unstable differential equation was treated: 


(27) —=—; y= 


tol 


at %4= 


vl 


to be integrated from x = 1/2 tox = 1 with hy = 2“ and 6° = 2°”. Results are 
given in Table 4. 

This illustrates the quality of the starting process in keeping the early truncation 
error small, a very important consideration in this case because such early errors 
are ultimately magnified one millionfold. Six significant decimals are preserved 
correct through 63 steps, in each of which the solution increases by 25 per cent on 
the average. The final error exceeds 2” because significance cannot increase. 

Each of the above tests required only 3 to 15 seconds of computer time, and 
some sort of longer test seemed appropriate. As such the author chose Bessel’s 
differential equation of order 16, and in particular, to find Ji(z) by integrating 
from z = 6 to z = 6000. In this range the function begins very small, increases 
monotonically and rapidly over 200,000-fold, and then makes almost 1000 complete 
oscillations. We put z = 2"z and hy = 2“ and 8° = 2°” and 2™ respectively, for 
two tests. Tables 5 and 6 show the results of this computation. 








TABLE 4 
Steps | x y Error 
0 .50 .000 000 476 837 
1 .507 812 5 .000 000 650 187 .O"1 


63 1.0 .500 000 546 694 .000 000 55 
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TABLE 5 
(6 = 2) 
Steps z J (2) Error J’ (z) | Error 
re 
0 6.0 | .0°1 201 950 | 0% 986 480 | 
1 6.125) .0°1 633 713 | .O41 | .0°3 963 765) .01 
77 | 16.0 | .177 453 370| .0% 177| .062 487955 .0% 066 
93 18.0 .261 082 210 | .05% 266| .003 519 524 .0% 005 
109 20.0 | .145 179 990 | .0% 150 |—.116 956 059 | —.09 118 
49 005 | 6132.0 | .004 126 972 | —.053 498| .009 311 583 | 
49 021 | 6134.0 | .006 748 858 | —.05O 809 |—.007 627 018 


49 037 | 6136.0 |—.009 741 657 | .0°4 174 |— .002 961 758 
49 053 | 6138.0 | .001 359 819 | — .0°2 666 | .010 089 144 























TABLE 6 
(G~ = 2%) 
Steps | z J (2) | Error | J’ (2) Error 
0} 6.0 | .95 201 950 | | .0°2 986 480 | 
1 | 6.125] .051 633 713 | 0" 03 963 765 | .0"1 
99 | 16.0 | .177 453 297| .0% 104| .062 487 925 | .0% 036 
126 | 18.0 | 261 082 096 | 0% 152| .003 519 520) .0% 001 
156 | 20.0 | .145 179 923| .050 083 |—.116 956 010 | —.0%O 069 
98 709 | 6132.0 | .004 130 418 | —.0% 052| .009 314 0689 | 
98 741 | 6134.0 | .006 749 685 | .05 018 |—.007 631 186 
98 773 | 6136.0 |—.009 745 792 | .0% 039 |— .002 960 774 | 
98 805 | 6138.0 | .001 362 434 |—.0% 051 | .010 092 495 








Some of the properties of the automatic interval control are well illustrated by 
these two tables. In spite of our asking for less than full register accuracy, the com- 
puter starts accurately enough and with a small enough interval in both cases so 
that the initial truncation error is half the least count, for it recognizes via test 
(22a) that early errors may be magnified by the instability of the differential equa- 
tion itself. The ultimate error is somewhat but not much larger than asked for, as 
it must be expected to be because of significance considerations. The interval is 
halved over most of the range and the error drops by just about 2~° as between 
Table 5 and Table 6 (due allowance being made for the change of phase of the 
error between the two calculations). In the calculation of Table 6 we end up with 
almost as many correct significant figures as were given initially. A further increase 
in e (and in computing time) would presumably improve the preservation of signifi- 
cance a little more. 

We emphasize that the above treatment of the Bessel equation is not claimed to 
be a good way of calculating Bessel functions, but was chosen purposely to illustrate 
how the method handles a rather “‘ill-behaved”’ problem. 
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Appendix A. The working formulas and truncation errors for degrees m = 2 
through 6 are collected here. 


m = 
y(z + h) = y(x) + hff(z) + a(x) + & [f(a + A) — fr)} 
(Al) f= f(z) + 2a(z) 
a(t +h) = a(z) + $[f(z + h) — f?] 
h* 
E, ae 1 4 ii yi" 
m= 8 
y(z + h) = y(z) + hif(z) + ale) + d(x) + 3 Uf(e +h) — fr)} 
f= (x) + 2a(z) + 3b(z) 
(A2) a(z+h) = a(x) + 3b(z) + 2 [f(x + h) — fr] 
b(ix +h) = b(z) + afi + h) — fr) 
~~ wa 
E, = (25/6) = y" 
m= 4 
y(z +h) = y(x) + hif(e) + az) + d(x) + c(x) + F8hLf(e + h) — fr} 
f= f(z) + 2a(xz) + 3b(x) + 4c(z) 
(a3) a(z +h) = a(z) + 3b(z) + 6c(x) + LLf(e + h) — fr] 
b(t +h) = b(xz) + 4c(x) + 3[f(x + h) — fr 
ee +h) = c(z) + sulf(z + h) — fr) 
6 
B, = (27/2) % y™ 
m=6 
y(z + h) = y(z) + hif@) + a(z) + bz) + c(z) + d(x) + Pelf(e + h) — fr)} 
f= f(x) + 2a(x) + 3b(xz) + 4c(xz) + = 5d(z) 
a(jz +h) = a(x) + 3b(x) + 6e(x) + 10d(xr) + $4[f(x + h) — fr] 
(A4) b(2 +h) = b(z) + 4c(x) + 10d(x) + S8[f(e + h) — f?] 
e(z +h) = e(x) + 5d(x) + #elf(z + h) — fr] 
d(z +h) = d(x) + sholf(z + h) — fr) 
E, = (863/12) J yv 
G3 
m=6 


y(z + h) = y(z) 
+ h{f(z) + a(z) + bz) + c(z) + d(x) + = e(z) 
+ s8285Lf(@ + h) — fr)} 
f= J(z) + 2a(xz) + 3b(x) + 4c(z) + 5d(x) + 6e(z) 
a(z + h) a(z) + 3b(x) + 6e(x) + 10d(xz) + 15e(z) 
+ 433i + h) — fr) 
b(x) + 4c(x) + 10d(x) + We(zx) 
+ SLi + h) — fr) 
c(z) + 5d(x) + 15e(z) 
+ Sif + h) — fr) 
d(jz +h) = d(z) + 6e(z) 
+ golf(z + h) — fr) 


(A5) b(z + h) 


c(z + h) 





an 


1} 


’ 
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e(x +h) = e(z) 


+rho e+ h) — pr 
inn 
E, = 5135 y"™" 


Appendix B. The flow chart (Figure 3) presented here is probably in terms which 
are general enough to apply to most stored-program computers. As shown, it pro- 


initial entry reversing norma! 
(supply x_, y°) entry entry 


gh: + 


e 
reset h,; construct 8 “/h | delay doubling h? 







































































































































































































°. 
save xX» Y¥i3 note n ¥ yes 
c 7 v | reduce delay 
reset stepping switch; 
clear a, b, c, d | tests oversatisfied? 7 
- yes 
| auxiliary subroutine 7” 41 ' [ [hy =h.? } >> 
4 r + 
0 
reverse h; 4- | reverse h_ }>—@ _¥ 
step doubling > | & - x.) = 0 (mod|2h] )? | be 
delay — 4 ¥ yes 
ot | double h ny 
A a 4 Lpl advance x and predict bt 
we s ¥ 
t | auxiliary subroutine 
ia e 
= ia 7. 7: —_ 7 
| iterate implicit equations | 
° tt P 
| restore y; 8 i gre— once 
we n 
g | record test results | A 
pm Spon 
Ls a s se 2-24 
RI J ot 12 ' did test (22a) rmiE 
fail t ’ 
fail gees 1-2 
| st¢ test (22b) sf ¢ 
fail? a h [did test (22b) fail? 
i wfail 
halve h; stop avs, - oy 
10's cate tiees 16 de-advance x; de-predict | 
At — h; stop if h 
under flows 
Aa _— 
nfl 29 |correct a, b, c, d 
az 4 
“= 
Pr guard 
digits =1/2 pe 





Idouble h ie 
(run) 


exit 
Fic. 3.—Flow Chart for one elementary step of integration. 
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vides for one elementary step of integration per entry into the routine, so that a 
master program can supervise the general course of the computation with complete 
flexibility. It also appeals to an ‘“‘auxiliary subroutine” (closed) to calculate f(z, y) 
given z and y, for complete flexibility as to what system of differential equations is 
being treated. The parameters which must be supplied are: the order n, the location 
of the auxiliary subroutine ho , the accuracy parameter e, and the location of a 
working storage of 2 + 10n memory locations. The working storage contains a 
location for zo , one for x and, for each i(i = 1, 2, --- n), 10 locations containing 
respectively y; , fi, ai, 6; ,c; ,d; , guard digits for y; , f;”, y,”, and y; at the beginning 
of the current step. The location normally containing guard digits is used for pre- 
serving the initial y,’ during starting. At the conclusion of the starting process this 
location is set to } so that when a double precision increment-addition is made to 
y; , the normally rounded y; will appear in the first of the 10 registers for the use of 
the auxiliary subroutine. 

The computing time per normal elementary step in this method is about 30n 
multiplication times (21n milliseconds on the Illiac) plus twice the time required to 
calculate the derivatives. There are 6n actual multiplications performed, the re- 
mainder of the 30n being accounted for by additions and ‘“‘housekeeping”’. Abortive 
integration steps, i.e. those partially done and then undone because of test failures, 
require only 2n actual multiplications but about 20n multiplication times plus 
twice the derivative calculation time. The starting process is clearly the equivalent 
in time consumed of not less than 24 normal step times. 

These figures are for a computer without special address modification features, 
and the housekeeping time may be expected to be rather less where address modifi- 
cation features are available. 


Appendix C. Here we discuss the choice of rounding procedures to guarantee 
against persistent noise induced by rounding, in an otherwise stable iterative arith- 
metic process, i.e. a process producing a convergent sequence when applied in the 
real number domain. Although we are not able to state a general recipe guaranteed 
to work in all cases, we can cite a qualitative principle which clearly always tends to 
improve the persistent noise behavior and which leads to a guaranteed solution in 
our particular problem of rounding the multiplications in (A4). 

If multiplications and divisions are rounded in the normal way, namely, by 
replacing any number which is a fraction in terms of the least count, by the nearest 
integer in terms of the least count, we do not in general get the resulting sequence of 
integers converging to a unique limit, as can be seen in terms of a simple example. 
Consider the process x,4; = maz, + b, where the z’s are real numbers and m and b 
are constants with | m | < 1. The sequence {z,} obviously converges and converges 
to b/(1 — m). When iterative processes of this sort are done numerically the limit- 
ing value is not generally known in advance, the objective of the process being in 
fact usually to find the limiting value. Accordingly we reformulate the problem in 
such a way that b does not appear: let y; = 2, — tn; , 80 that y, obeys yn = MYn 
and tends to zero. Then we further reformulate so that the quantity eventually to 
be rounded in some sense, is the change in magnitude of y in an iteration. Spe- 
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cifically, we write 
(C1) Yori = £Yn F (1 — | m/)yn = +Yn F wYn 


according as m > 0 or m < 0. Observe that 0 < » < 1. 

The digital (integer) process corresponding to the real number process (C1) 
involves rounding the product yy, to an integer according to some rule. Using an 
asterisk to denote a quantity integral in terms of the least count, we have for the 
digital process: 


(C2) yn = £(yn* — [yn*)) 


where [ ] means some sort of rounding. 

Normal rounding causes most of the sequences generated by (C2) to misbehave. 
If » = } — e for example, then it is easy to verify that under normal rounding 
rules an initial yo* = 0 leads to the sequence 0, 0,0, --- ; initial yo* = 1 leads to 1, 
+1, 1, +1, --- ; all other positive initial yo* lead to 2, +2, 2, +2 --- ; and simi- 
larly for negative initial y*. This general sort of misbehavior is not peculiar to the 
value of u chosen for illustration, but is typical of most u’s. In the formulation (C2), 
however, the source of the difficulty is easy to discern: it is merely that the term 
{uy,*] normally rounded may often vanish when y,* does not, so that the magnitude 
of yn* may “get stuck” at a non-zero value. 

The difficulty is entirely removed in this simple example by redefining the 
rounding process so that 


P for x exactly integral 
(C3) [x] = {integer nearest (x + 4) for z positive non-integral 
integer nearest (x — 4) for x negative non-integral 


We term this special kind of rounding “rounding away from zero,” for it consists 
of moving the number z away from the origin just far enough to make it integral. 
So defined, [uy,*] does not exceed y,* in magnitude, is of the same sign as y,* and 
does not vanish unless y,* vanishes. Thus, all integer sequences generated by 
(C2) must now converge to 0. 

The general principle is accordingly that if we can formulate an iterative digital 
process so that the quantity to be rounded is a correction subtracted from the 
previous value of an integer variable intended to converge to zero, as in (C2), then 
the quantity to be rounded should be rounded generally away from zero. In more 
complicated cases where several integer variables are involved the correction (in 
the above sense) to each may be a function of all the variables; but still it should be 
rounded away from zero. 

Our particular problem consists of rounding the multiplications A[ |], B[ |, 
C[ ],D[ | in the working equations (A4). Suppose that f tends asymptotically to 
a constant and consider what may happen when a --- d have become small. 
Then f(z + h) — f(x) will cancel out of (A4) at some stage, and thereafter the rel- 
evant equations of the process will be: 
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6c43 = au? + 3b,,* + 6c,* + 10d,,* 
— [(25/24)(2a,* + 3b,* + 4c,* + 5d,*)] 


* 


baa = ° + 4c,* + 10d,,* 
— [ 35/72 (2a,* + 3b,* + 4c,* + 5d,*)] 
(C4) , 
Can = o,* + 5d,* 
—[5/48 (2a,* + 3b,* + 4c,* + 5d,*)] 
dnt sual d,* 


— [ 1/120 (2a,* + 3b,* + 4e,* + 5d,*)] 


where the asterisk signifies a quantity integral in terms of the least count, and the 
[ ] symbolizes rounding. Note that these equations are in just the form we require 
to apply the ‘‘rounding away from zero” principle, since the terms 3b,*, 4c,* etc. 
are integral and have no effect on the behavior of the rounding. 

Normal rounding in equations (C4) leads to persistent roundoff noise. The 
rounding process is so non-linear that we have no analytical theory and must work 
out specific numerical examples. Two examples of indefinitely persisting (cyclic) 
roundoff noise are: 





n a’ b* c a* n a* b* e da* 
0 1 0 0 0 0 0 1 0 1 
1 —1 —1 0 0 1 5 7 4 l 
2 1 1 l 0 2 6 8 4 1 
3 1 0 0 3 5 6 3 1 
4 —1 — —] 0 4 4 6 3 1 
5 —1 _ 0 0 5 5 7 4 l 
6 1 1 0 6 6 8 4 l 
ete. ete. 


As we saw in Section 8, any behavior like this (and there are many cases of it) can 
frustrate the interval control in its attempts to increase the interval when the 
interval obviously ought to be increased. Curiously enough, the persistent cycles 
of roundoff noise contribute practically no error to y, for the contribution to y, 
averaged over a repetitive noise cycle, is no more than about h/60 times the least 
count. However, proper behavior of the interval control alone is enough reason for 
rectifying the roundoff behavior. 

The simplest change in rounding which suggests itself is rounding all four multi- 
plications in (C4) away from zero. However, such a simple remedy does not work, 
for it represents too drastic a modification of the fourth equation of (C4). It im- 
plies in fact that d* must change unless (2a* + 3b* + 4c* + 5d*) is zero, and per- 
sistent oscillation of d* results inevitably. After some experimentation the author 
has concluded that the best rule is: round the first three multiplications in (C4) 
away from zero according to (C3), but for the fourth multiplication move the 
multiplicand (2a* + 3b* + 4c* + 5d*) away from zero by 16 unitsand then multiply 
by y4s, rounding normally. The treatment of the fourth multiplication is a “partial” 
rounding away from zero or a less drastic modification of normal rounding, but 
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clearly in the same spirit. The rounding rules thus finally fixed upon will cause 
every initial quadruple of integers to converge to (0, 0, 0, 0), as was verified by 
letting the computer treat every case. Actually, all initial quadruples of integers 
between —2 and 2 inclusive were examined, and all tend to (0,0, 0,0). The average 
number of steps to arrive at (0,0, 0,0) is42 and the maximum is 14. If we move the 
multiplicand of the last multiplication only 12 units instead of 16, one persistent 
cycle appears. If we move it 14, 16, respectively 18 units all quadruples converge 
to (0,0, 0,0) but the average number of steps to clear begins to increase. Thus 16 
seems a safe compromise. 

These principles may be of help in deciding how to round the arithmetic in 
other iterative digital processes, such as solving systems of implicit equations. In 
our present state of knowledge of the subject a certain amount of experimenting of 
the sort described above will probably have to be done in every individual case 
more complicated than the one-variable case. The general reason for stabilizing 
roundoff noise in these ways is to improve the functioning of tests-for-end, for such 
tests are subject to the same difficulties as test (22b) in our procedure. 


University of Illinois 
Urbana, Illinois 
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Triple Product Integrals of Laguerre Functions 


By J. Gillis and M. Shimshoni 


1. Introduction. We shall use the following standard definitions for Laguerre 
polynomials (1) and Laguerre functions (2): 

(1) Lx(z) = 2 (—1)' (*)2 
r=0 ry Tri 
(2) Aa(z) = € *"L,(z) 

The Laguerre functions are known to constitute a complete orthonormal set in 
L?(0, ~ ). Given a differential equation over 0 < x < ~ one naturally thinks, there- 
fore, of the possibility of solution by expansion as a series of Laguerre functions. 
However, for this to be useful for non-linear differential equations, we need to be 
able to expand the product of two Laguerre functions as a linear series of these 
functions. The main part of this paper is devoted to methods for effecting the 
expansion, and we shall also give an application of the results. The similar problem 
for Laguerre polynomials has been solved by Watson [5], and methods for com- 


puting the expansion coefficients discussed by Gillis and Weiss [3]. 
We may write 


(3) ArAs = > CrstXt 
t=0 
where 
(4) Cas = [ Ar(a)Ag(a)Ag(x) dz. 
0 


We shall also write this as C,.,.;. 

It follows that the coefficients will be symmetric in all three suffixes. We give 
below a table of these coefficients for 0 < r S s < ¢ S 10, and also expressions of 
C,s: as polynomials int for0 S r Ss S 3. 

In Section 2 we discuss a number of formulas for C,.,. Those in (a) and (b) 
involve three-fold summations and apply to general r, s, ¢. In (c) we obtain a 
simple sum formula valid for the case r = 0 and, in (d), a double sum for the case 
r = s. In Section 3 we shali derive two recurrence formulas for the coefficients 
Cs: . As a check on the stability of these latter formulas in practice, it is advisable 
to have comparatively simple alternative methods for computing C,,., for certain 
particular triads r, s, t. For this purpose the formulas of Section 2 (c), (d) can 
be of use. 


2. Explicit Formulas. 
(a) 


) r 2 Bt Y 
pe —3r/2 _1ye(7\2. _ 16 (8\ 2 (—1)7 ‘= 
Cont I eet (-1) (4x 1) (o> 1) (f = dz 


ys ie (a) (:) (;) (‘) (a+ 6 + 7)! 
3 a,B.7 3 a} \B/ \v alBly! 


Received May 8, 1961. This research has been supported by the Air Research and De- 
velopment Command of the U. 8. Air Force (European Office). 
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(5) 











le- 
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(b) Since 


x" = nd (—1)* (") L,(z), (ef. {1}), 


a=0 


it follows that 


Lictiis) p> (ar? "= yi *(°)2 


(6) 
p+¢@ 
a7 (= wT Va (-1'(? +9) nx), 
p=0 q=0 Pp q P —_ ‘ 
i.e., 
(7) L,(z)L,(x) = 3 Aretl( x) 
where 
a J + q P + q 
8 Ant = (—1)? «(’) (°) (? )( ). 
(8) (=22 p)\a)\ p 
Now [1], 
(9) as 2 = > ¥ 3° "L,(z) 
p=0 
and so 
4 Creel (2) = ie *L,(x)L,(2x) = > 8 . ‘A remlip( X) Lm (2) 
7 =0 m=0 
(10) z 2 m+p 
= 2>> bs } te ee mptl(2). 
p=0 m=—0 t=|m—p 
Hence 
2203 hay ‘A - A mpt = 2 > [nf Oe 
vail? rs m,p,a,B, 7,6 


MEME rs MNS TD) 
a)\8 a m v/\i Y t 
Pp q-r-1 _ oF! 
(?) gt = 
while the sum over m is 


mn (at+B\(/m\ _ ,_,\f«+8 4)" chili 
2 (1) ( m Mr) =‘ »( 7 NE ™ ( m— 


=(-1)” if y=at@B. 


The sum over p yields 


M Ls 


(12) 


Pp 


IV 


(13) 


Hence 


; “ae ae: (—1)stH48 r\/s\fat+6\f/at+64+6\/at+6+s5 
(14) Cue = 2 2 2 — — (oy @ i y \ ' ) 
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(c) Take Laplace transforms of both sides of (3). This gives [2] 


' ‘) p '(p—1) *" Fil—r, —s; —r —8;1 — p’] 
(15) 


; 
= 2 Culp — )'(p + 4). 
Writing g = (p — 3)/(p + 4), we obtain 
Tema =2(7F Jat ge — a 


; 
; eee et D 
al ‘eal | 


In the special case r = 0 this leads to the simple result that 


(16) 





Cost = coefficient of g‘ in 2(1 + q)*(3 — q) 


(17) = 23" "(6 + n) in! (t — n)! (8 +n — t)y7 
n=0 


We could similarly use (16) to obtain expressions for C,,, for small nonzero values 
of r, but the formulas soon become prohibitively complicated. It follows incidentally 


from (17) that Co, > 0 for all s, t and that 3°*'*'C,,, is an even integer. 
(d) From the generating function [1] 


(18) ¥ La(z)u" = (1 — uw) exp {ue(u — 1)4 

we obtain 

(19) Do da(z)u" = (1 — u) exp {He(u + 1)(u — 1) 
Hence 


(0+) p(0*) 
r-(2)Ae(2) = (2x) | / 6 eye + 9) 
-exp {3 | = + ~ z |} du dv. 


(21) (w + 1)/(w — 1) = (u+ 1)/(u — 1) + (0 + 1)/(v — 1) 


and write 


(20) 





Choose w so that 


exp {3(w + 1)(w —1)7} = (1 -— wo) diz) 
Then 


x (0+) p0(*) 
(22) d-(z)ae(2) = 32 (Qui) rj (z) / | F, (us, 9)u710-** dude 
t=0 











TRIPLE PRODUCT INTEGRALS OF LAGUERRE FUNCTIONS 53 


where 


(23) Fi(u,v) = 2.1 +u+v — 3w)‘(3 — u — vo — w) 


It follows that C,,; is the coefficient of u’v’ in F,(u, v). 
By putting u = 0 we can calculate C,,, and immediately recover (17). Formulas 


r 





0 
for other small values of r can be obtained by expanding [ au Fu, ®)| in powers 
u=@ 


of v. However, there is another case in which this can be used to obtain a closed 
expression for C,,; , namely, when r = s. For this purpose we write 


(24) Fu, v) = (—1)'2* "(4a — B)*(1 — 38)" 
where 
a= w 
and 
(25) B = (u+1)(v + 1). 


We note that the coefficient of u’v’ in the expansion of ” is ("). But, by (24), 


(26) Fu(ae) = (—1)'F** YY (-1)"(F)i4aymarrae(*t "). 


m=0 n=0 
It follows that C,,,, the coefficient of u’v’ in the expression F,(u, v), is 
t CJ 2 
97’ YY 7 t __1)™g-2im+tn) t t+n m+n 
(27) 4(—1) 2, 2, | ‘ys (*)( n tee F 


The table of numerical values of C,,, suggests the conjecture that C,,, has the sign 
of (—1)‘ for 0 < t < r. This would mean that the double sum in (27) is always 
positive for 0 < ¢ < r, but we have not been able to prove the conjecture. 


3. Recurrence Relations. 
(a) A Two-Index Recurrence Relation. From the relation 


(28) rrA(z) = —(r + 1)Arga(z) + (2r + 1)A( 2) — rAu(2), 
[1], we deduce that 
(29) [ te(22)desa(2) dx = —(r + 1) 
and 
(30) [ rh, dx = 2r +1, 
0 
while 


(31) [ ad, = 0 este.r dt 
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Multiplying both sides of (28) by A, and using (3), we get 


Zz. CrettXt 
t=0 


= —(r + 1d Cr41,8,0X4 + (2r + 1) Do Crd = r 2s Coase 


We now multiply both sides of (32) by \, and integrate from 0 to ~, using (29), 
(30), (31), and immediately obtain 


(r + 1)Cr4i.8,t 
(33) 
- ebcs + 2(r — C4. + (t + 1)C,.2,044 pos eee - 
This is the required recurrence relation. It can be used quite effectively, in con- 


junction with (17), to compute a table of C,,,. We first compute Co, by (17) for 
an adequate range of (s, ¢). It then follows from (33) that 


(34) Cist = tCo,*,11 — 2tCosr + (t + 1)Co,s,e41 


and this gives us C;,; , etc. The computation is reasonably stable, although, as a 
safety precaution, one would carry more digits than are actually needed. Actually, 
the explicit formula for C,,, , obtainable from (17) and (34), is 


t 


(35) Cw: = 2-3" *(s — 14+ n) init — n)!(s —t+n+ 1)3 "Pree 
n=0 
where 
Fuse = (48 + 1)n* + (88° — 12st + 68 + 1)n 
36) 3 
( + s[4s° — (12t — 5)s + 9t(t — 1) + 1). 


A similar, but much more complicated formula can be derived from this for C2,; . 
The corresponding formulas for C,,, become very complicated for larger values 
of r. 

For work with an electronic computer it would be better to have a method for 
generating the C,,, as required rather than to store a table of these coefficients. 
In the next paragraph we propose a method of achieving this by a recurrence 
relation which operates on only one of the indices. 

(b) An Alternative Recurrence Relation. We recall that , satisfies the differ- 
ential equation 


(37) tre” + r- + (r + 3 — je)A, = 0. 
Hence, if we write 

(38) ur=2'r,, 

we can easily verify that 


(39) Ue” + Pit = 0, 
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where 





(40) p, = 4.4S%8.! 


r 


47? x 4° 


It may then be deduced, by a standard procedure, that y = u,u, satisfies the differ- 
ential equation 


ay? + a’y” + 2°(1 + 2ox — 2*)y” 


(41) 2 
—a(2+ ox +2°)y¥ + (2+ or+ &2)y =0 

where 

‘iis o=r+s8s+1) 


6=|r-—s| j. 

Let z = AA, . Then, by (38), 

(43) y = zz 

and this may be substituted directly into (41), giving 

(44) 2&£(z) = a2'2°" + 5az” + (4 + Qor — 2°) 2” + Blo — x)’ + (8 — 1)z = 0. 


We now write 
(45) z= p> CretXt 
t 


and substitute in £(z). 
After a little manipulation, using the properties of Laguerre functions, we obtain 


(46) £(X) = (« - p — 2) n +[¢ — 2p +8 -5 475" r—3a'|n 
where 
=t+}3 
and hence, again going back to basic properties of A, , 
rL(AL) = Pelt + 1)(t + 2)(E + 3)Aras 
— pe(t + 1)(t + 2)(10t — 86 + 13)Arse 
—s(t + 1)[30 + 2t + 2(88° + 40 + 1))rcus 
+ (§(2t+ 1I)(@ +t+1) + H+ 1) — of +#4+ 4) 
— ppt(30 + 4t + 168° — 8o + 3)A4 


(47) 


— pet(t — 1)(10t — 80 — 3)Ave + Pet(t — 1)(t — 2)AL3. 
It follows from (45) and (47) that 


(48) t&£(z) = > Ba. 
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where 
B, = x(t + 3)(t + 2)(t + 1)Cr 0.145 

— pe(t + 2)(t + 1)(10t — 86 + 17)Cy 0142 
— ze(t + 1)[30 + 10 + 168° — 80 + 10]Cy 6.141 

(49) + (§(2t+ 1)(@ +441) + P(2t4+ 1) — off +t4+ 4))Cu 
— pet[3e? + 4t + 168° + 80 + 3)Cr e104 
— pet(t — 1)(10t — 80 — 7)Cy2,1-2 
+ Fot(t — 1)(t — 2)Cr 0-2. 

Since £(z) = 0 we have, from the completeness of the \,’s, that 


C _ 10t — 86 + 17, 4. 38 + 10t + 163° — 80 + 10, 
47,8 t+3 — 3(t + 3) /T 8 t+2 3(t re 2)(t + 3) ~——— Vere ti 
_ 10(2# + 3@ + 3t + 1) + 26°(2t + 1) — 160(¢ + t + 3) 

3(t + 1) + 2)(¢ + 3) 








~) 
Cea 








(50) 2 2 
4 U3t — 4¢+ 166 + 8¢ + 3) 
3¢+ Dt+2¢+3) 
4 Ki — 1)(10¢ — 80 — 7) i iit — 1)(t — 2) ’ 





~~ ; Cra, —3- 
3(t + 1) + 2 + 3) (+ 1) + 204 3) 
Suppose, then, we know C,,, for ¢ = 0, 1, 2. Then, putting ¢ = 0 in (50) gives us 


18C,..3 = —2(8¢ — 17)C,,22 
(51) . 

+ 2(88 — 40 + 5)Cr21 — (28 — 86 — 10)C,...0 
and the subsequent terms are all obtainable from the recurrence relation. C,,¢,0 
and C,,,,, can be computed directly from (17) and (35) respectively. It would be 
possible to develop a corresponding formula for C,,, » but it would not be very useful. 
Perhaps the best way to obtain C,,,.2 is by use of the relation (33). 

It should be remarked that equation (50) is probably the most suitable, among 
the formulas given above, for use with an electronic computer. For work with a 
desk computing machine it is rather complicated, and there is little doubt that (33) 
would prove more useful. 


4. Tables. When working by hand, it will generally be convenient to have a 
table of the numerical values of the C,,, . We give this immediately below in Table 
Ifor0 < r= s St S 10. Some general formulas for C,., as polynomials in ¢ for 
0 Srss S 8are given in Table II. 


5. Application. As mentioned in Section 1, our purpose was to apply the above 
ideas to the solution of non-linear differential equations over a semi-infinite range. 
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TaBLeE I 
| | i | 
a s t Crat r ne t Crst 

0 0 0 |  — .66666667 0 | 5 | 8 |  .07837182 
0 0 1 | 22222992 Oct | Bal .04878490 
0 0 2 | .07407407 0 | 5 | 10 | .02839004 

0 0 3 | .02469136 | 
0 0 4 | .00823045 0 6 | 6 | . 16045055 
0 0 5 | .00274348 0 ) . 14305131 
0 0 6 00091449 0 ei @ol . 11236326 
0 0 a .00030483 0 6 | 9 | .07984000 
0 0 8 .00010161 0 6 | 10 | 05233831 

0 0 y | .00003387 | 
0 0 | 10 .00001 129 0 £1 2x1 . 14885106 
0 7] 8 | . 13475521 
0 1 1 | . 37037037 0 ? | Syl . 10898616 
0 1 2 | . 22222222 0 7 | 10 | .08038815 

0 1 3 . 10699588 | | 
0 1 4 .04663923 0 8 8 . 13945383 
0 1 5 | .01920439 0 8 9 | .12773173 
0 1 6 | .00762079 0 8 | 10 | . 10570201 

0 1 A .00294671 | 
0 1 : | 00111772 0 9/ 9] . 13163910 
0 1 9 .00041773 0 9 | 10 | . 12169095 

0 1 10 | .00015430 | | 
| 0 | 10 | 10 | . 12500700 

0 2 2 . 27160494 | 
0 2 3 |  — .20027435 1 | 1 | 1 | —.07407407 
0 2 4 | . 11796982 Bese . 17283951 
0 2 5 | 06127115 1 | 11] 3 |  .21124829 
0 2 6 | .02936544 ick | @o . 15089163 
0 2 . | 01331098 Bigs cl 5 | .08687700 
0 2 8 | .00579180 1 | 1] 6 |  .04440380 
0 2 9 | 00244242 bok Po .02103338 
0 2 | 10 .00100482 1} 1] 8 | .00944978 
1/1] 9] .00408324 
0 3 3 .22405121 edcod] wal .00171233 

0 3 4 |  .18076513 | | | 
0 3 5 | . 12000203 1 | 2] 2 | —.04115226 
0 3 | 6 | 07021287 | 1 | 2 | 3 | .08504801 
0 3 a .03762977 Thus i £2 . 16369456 
0 3 | 8 | — .01891085 1 | 2| 5 | [15333029 
0 3 9 | .00904835 cle | @al . 10841843 
0 3 | 10 | .00416520 Soke 1 2 06553879 
1 | 2 8 .03580078 
0 4 4 | . 19519382 Beading 9 .01821086 
0 4 5 | 16532033 1 2 | 10 .00878492 

0 4 6 | .11851174 | 
0 4 Tol .07545146 .. 3 | —.02042372 
0 4 8 | .04399736 Birdy. @ 4 | .05009399 
0 4 9 | 02398552 Ricoh rt 5 | . 12508256 
0 4 | 10 | .01239969 Bok if 6 . 14086606 
th tn 7 . 11596019 
0 5 5 | . 17527816 Te ear . .07989974 
0 5 6 | . 15303674 et a 9 .04897920 
0 5 7 | 11566665 1 | 3 | 10 .02762203 
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TABLE I—Continued 




















| s t Crst r $s t Crest 
1 | 2 a. — .01188843 Spy 10 | .07467760 
1) 4 5 | .03383631 | 
1 4 6 | .09794126 ee 4 | . 12284713 
1 = 1 if . 12524062 ee 5 | .04609731 
1 =) 2 . 11552364 psig 6 | —.00911108 
Pickle | 9 .08823376 | 2 | 4 | 7 | — .01571953 
1 4 | 10 .05941345 Bf 8 | 06735022 
| | BO 9 | .09906079 
1 5 | 5 | —.00795949 ek ee 10 | . 10112375 
1 5 | 6 | .02489458 
1 S|. ie .07889869 2 | 5 5 . 10347339 
1 5 | 8 | .1104813 | 2 | 5] 6 .05127568 
1 5 | 9 | . 11111857 2 | 5 7 — .00230192 
I 5 | 10 .09212795 2) 5 8 .00330548 
2 | 5 9 04431534 
1 6 | 6 | —.00582692 Sg | .08081952 
1 ABA .01934154 | 
1 6 | 8 | 06515827 
1 6 | 9 | .09753256 2 6 6 | .09128849 
1 6 | 10 | . 10499757 2 6 7 | .05266394 
2 6 | .00397312 
1 7 | 7 | —.00451393 2 6 9 | —.00198035 
1 7 | 8 | .01560049 | 2 6 | 10 02815229 
1 7 9 | .05492763 | 
1 7 10 | 08655013 2 7 7 08275543 
| 2 7 8 | .05253757 
1 8 | 8 — .00363396 2 7 9 .00903873 
1 - | 8°] .01293342 | 2 7 10 — .00363722 
1 8 | 10 04709154 
2 8 8 .07631313 
Sei | 9 — .00300868 2 8 9 | .05177168 
respec | 10 | 01095104 2 8 10 | .01298846 
1 | 10 | 10 | —.00254500 | 2 9 9 .07120551 
| | | | 2 9 | 10 | .05073310 
nis at 2 20576132 | 2 10 10 .06701828 
2 2) 3 — .00457247 | | | 
2 am .00335315 a 1 3 3 — .05375197 
Scie |. 5 .09053498 an 4 . 11295704 
Sree | 6 . 13260174 3 3 5 07694175 
Pog 7 . 12277939 3 3 6 — .00123814 
puree | 8" .09060272 3 3 7 .00106252 
2 B ind] .05811997 3} 3 8 .05014082 
2 oie .03385805 Be ryiry 9 .09001406 
| PVigtoeetog 4 . 10048996 
2 Kase . 15658182 
2 R fond .03119443 ie 4 — .03744912 
2 Pg — .01249810 3 | 4 5 .07177090 
2 3 6 | 04167161 3 | 4 6 .08787554 
2 3 | 7 | 09805416 3 4 7 02244924 
2 prope ng ey . 11549354 Seerig 8 — .00731554 
Sg | 8°) -10130123 | 3 | 4 9 


.01802366 


| 
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TaBLeE I—Continued 
r | s t Crat r s t Crat 
| | — + =e ena 
3 4 | 10 .05982871 4 7 10 .02567632 
3 5 | 5 .02226023 4 8 8 .06531960 
3 5 | 6 .04824409 4 8 9 .02633548 
3 5 | 7 .08353180 4 8 10 — .00693875 
3 5 | s 03995357 
3 5 | 9 00156912 4 9 9 .05946498 
3 | 5 | 10 .00178650 4 9 10 .02814148 
3 | 6 | 6 01411745 | 4 10 | 10 05504858 
Aw ee .03516063 
yt ee ig 07516477 5 5 5 — .04386690 
3 | 6 | 9 .04998029 5 5 6 .08650875 
3 | 6 | 10 .00808104 5 5 7 .03500443 
5 5 8 — .01294078 
Por res 00985944 5 5 9 .02870872 
Va Pee ee .02724675 5 5 10 .06403893 
axe 06668946 
- | tis 05474855 5 6 6 — .03350414 
| | 5 6 7 .06117194 
3 | 8 | 8 .00741466 5 6 8 05170139 
ie ee .02201080 5 6 9 — .00581844 
Tih, Go. e588 05915185 5 6 10 .00668642 
3 | 9] 9 .00586119 5 7 7 — .02189833 
3 | 9 | 10 .01830449 5 7 8 .04397838 
| | 5 7 9 .05640500 
3 10 | 10 00479475 5 7 10 .00564310 
4 4 | 4 . 12925988 5 8 8 — .01467602 
4.1344 4 .02134197 5 8 9 .03341818 
4 | 4| 6 .02196251 5 8 10 05551730 
7) <5 2 .08246301 
4 4 | 8 .05322789 5 9 9 — .01055417 
4 4 | 9 00474524 5 9 10 .02668086 
Qh ft 1 .00303408 | 
| | 5 10 10 — .00808528 
4 5 | 5 . 10673496 
4 5 | 6 00196740 | 6 6 6 .09687151 
go Pg hog .00100286 | 6 6 7 — 02410491 
Asti. nda of .06117005 6 6 8 .02789589 
eg, ieee x ee .06583263 | 6 6 9 06239762 
4 | 5 | 10 02450630 | 6 6 10 01401545 
a ST % .08660706 | 
2 ae ie ee .01580477 | 6 7 7 .08320655 
4 6 8 .00811483 | 6 7 8 — .00757894 
4 | 6 9 .04086770 6 7 9 .00739178 
4 | 6 | 2 .06525394 | 6 7 10 .05666489 
a| 5 eee @ 07367692 | 6 8 8 .06900813 
4 si. 02281235 | 6 8 9 .00401100 
4 7 9 } 6 8 10 — .00218959 


.00867086 
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TaBie I—Continued 



































r s t Creat r | s | i | Crst 
6 9 9 05911346 . | St a8 | .07864342 
6 9 | 10 01077477 8 8 | 9 | —.02407255 
ane | ot .02947714 
6 10 | 10 .05253399 | 
8 | 9 | 9 | — .06917670 
7 ek — .03776986 8 | 9 | 10 | —.01154288 
7 yi .07123199 
7 7 | 9 01652253 | 8 | 10 | 10 | 05829179 
7 7 | 10 — .00951384 
| 9 9 | 9 | —.03353927 
7 8 | . — .03035843 9 | 9 | 10 | 06114106 
7 8 9 .05351461 | | 
7 8 | 10 | — .03227809 ei ei toe 
| | 9 | 10 | 10 | —.02786041 
7 9 | 9 | —.02107272 
7 9 | 10 |  .04012854 | | 
| | 10 | 10 | 10 | = .06681991 
7 10 | 10 | —.01469083 | | | 
TaBLeE II 
r s Crst 
0/0] 2-3-4 
0 | 1 | (8t + 2)3-"? 
0 | 2| (162 + 2)3-" 
0 | 3 | (64 — 482 + 56t + 6)3-"5 
1 | 1 | (822 — 48¢ + 10)3-** 
1 | 2 | (64 — 2402 + 200¢ + 18)3-** 
1 | 3 | (256¢ — 1664 + 30562 — 1264¢ + 78)3-** 
2 | 2 | (128 — 1024 + 24642 — 1664¢ + 66)3-*5 
2| 3 | (512¢5 — 6528¢ + 28160 — 46848 + 24824 + 438)3-*7 
313 re 3993605 + 2858241 — 920832 + 1314560" — 647280t + 
4410)3-* 











An example of such an application will be given to illustrate the various technical 
problems involved. 


The Blasius Equation. 
(52) y” + yy” =0 
(53) y(0)=y¥(0)=0;  y¥(#) =2. 


This is a well-known equation whose numerical solution has long been known [4]. 
However, as an example of the Laguerre functions method, we sought an approxi- 
mate solution of the form 


(54) y = fy(z) = 22 + ay + ) 2 byrA, 


r=0 











al 
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where ay , byo, bvi,--- , bww were constants to be determined. 
The boundary condition at infinity is clearly satisfied, and, to attend to the 
conditions at z = 0, we need 


(55) ay + > be, =0 
and 
(56) 2- 2X (r + 4)bwr = 0 


thus leaving us the possibility of imposing N further conditions on the (N + 2) 
coefficients ay , by, (r = 0,1, --- , N). As an obvious set of conditions, we should 
substitute from (54) in (52), express y” + yy” as a linear sum of i,’s, and equate 
the coefficients of Ao, A1, --- , Aw-1 to zero. However, there is a difficulty in the 
way of this program. The expression for 4,’ in terms of the Laguerre functions 
themselves involves all the A, for 0 < r S n, while the expressions for \,.”, A,” 
become very complicated. The situation can be saved by pre-multiplying (52) by 
x and making repeated use of the relation 


(57) Try’ = 3(m + 1)Anqi — 9A, — 4mA,4 - 


Proceeding in this way, one finds oneself confronted with expressions of the form 
2h, . They, too, can be resolved by the repeated use of the relation 


(58) DrAn = —(n + 1)Angi + (22 + 1)A, — MAW, 


which also follows trivially from the fundamental properties of Laguerre poly- 
nomials. 

After making all of these substitutions into (52), we still have to deal with 
products of Laguerre functions arising from the nonlinear term, and these have to 
be resolved by (3). We are now in a position to equate the coefficients of A» , A, , --- , 
Aw_1 to zero, obtaining a set of N quadratic equations which, together with (55) 
and (56), should suffice to determine the coefficients. In general, there will be 
more than one solution, and any one of them might, if N is large enough, be ex- 
pected to yield a function fy(z) which will approximate the exact solution of (52). 
Since the solutions are obtained by equating the coefficients of \», A; , --- , Aw—a 
to zero, it has been found useful in practice to select the one for which the coefficient 
of Aw is least in absolute value. 

Setting up the equations even for so simple a case as (52) is not a trivial task, 
and can become extremely laborious for more complex equations. However, there 
would be no real difficulty in having all the work, including the formal steps repre- 
sented by (3), (57), and (58), programmed for an electronic computing machine. 

We have described the procedure so far in some detail, since it is of quite general 
application. The next step is actually to solve the equations for ay , by, (r = 0, 
1, --- , N), and for this purpose the following type of approach has been found to 
be practical. One first solves for some small vaiue of N. The advance from N to 
N + 1 is effected by solving the equations for N + 1 by the Newton-Raphson 
method, taking as a first approximation to this solution 


(59) ai: = av, O82, = bwr (r = 0,1,---, N), dS = 0. 


bwat = 








62 J. GILLIS AND M. SHIMSHONI 


An advantage of this choice of first approximation is that it might, for obvious 
reasons, be expected to be better as N increases. Hence the number of steps re- 
quired for convergence decreases. 

In the particular case of the Blasius equation, we started with N = 1, i-e., 
with three coefficients to be determined. Eliminating two of them by (55), (56) 
left us with a quadratic equation in one unknown. The step to N = 2 was effected 
as described. There would be no difficulty in principle in carrying on to higher 
values of N. However, the arithmetic soon becomes extremely laborious, and the 
task is best handed over to an electronic computing machine. The result for N = 2 
is shown in Table III. The function tabulated as f(z) was obtained by direct 
numerical integration, using essentially Hamel’s original method, and is given cor- 
rect to three decimal places. 


Taste III 

x Sa{x) f(z) 

0 0 0 

1 0.732 0.650 
2 2.337 2.305 
3 4.247 4.280 
4 6.252 6.279 
5 8.248 8.279 
6 10.273 10.279 
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The Behavior of Pseudo-Random Sequences 
Generated on Computers by the Multiplicative 
Congruential Method 


By V. D. Barnett 


1. Introduction. The use of pseudo-random elements in Monte Carlo work is 
essential when the scale of this work is such that the calculations involved are too 
extensive for hand calculating machines and it is necessary to employ an electronic 
computer. Although the ability of modern digital computers to perform simple 
binary operations at very high speed makes their use in this work particularly 
relevant, the limited extent of the computer memory, the relatively slow input 
speeds and speeds of access to the memory, and the very large number of random 
elements required (often of the order of 10°) combine to make it unfeasible to pre- 
pare the elements beforehand in the required form for input to the computer (e.., 
on tape or cards), and we must resort to some means of generation of the random ele- 
ments within the computer. 

Mechanical means of generation on peripheral equipment, e.g., by radioactive 
decay, thermal noise in electronic valves, etc., are undesirable because of the irre- 
producibility of the numbers obtained, which enables no check to be kept on their 
quality. It is therefore natural to employ some deterministic method of generation 
of the random elements by recurrence relationships. One such technique which has 
attracted much attention is the multiplicative congruential method (see for exam- 
ple [1], [2] and [3]) which proceeds as follows: 

Choose p, x at random as the starting values. 

Successive random elements are then obtained by the recurrence relationship 


(1) 2r41 = x-2z, (mod M) 


with zo = p: that is, 2,4; is of the form pz” (mod M). 

Using this method to generate pseudo-random sequences, one must place 
certain restrictions on x and p to ensure that the process does not degenerate to 
zero and that the maximum possible cycle of distinct elements is obtained. We 
shall see that only two restrictions must be placed on x and p to produce the maxi- 
mum cycle and that it is also possible to describe fully the behavior of the system 
for x and p not satisfying these restrictions. 

Because of the binary base of most digital computers, it is convenient to choose 
M of the form 2*. Reduction, modulo M, is then simply a shift of the product 
px’ to retain the least significant k binary digits. If the computer, in common with 
many modern computers, has a facility for low multiplication, i.e., multiplication 
which retains only the least significant half of the set of binary digits comprising 
the product, and the numbers are stored as ko digits, it is a further advantage to 
choose 


M = 2", 
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and this enables the successive random elements to be obtained by just one opera- 
tion. 

The length of the cycle of iterates and the conditions under which the maximum 
cycle is obtained must obviously depend on p, xz, and M = 2“, and although this 
system (for M = 2*) is discussed in the literature [1], [4], [5], there appears to be 
some confusion as to the behavior of the system for given p, z, and M. 

It has often been assumed desirable to choose x of the form 


2 2k+1 


z=5 
me AK +1 
z= 7* 
13 
x = 13° ete. 


to ensure the maximum cycle; see [3]. For instance, Taussky and Todd [1] state 
that for the system 


p =1 (mod 5) 
x is 5" 
M = 2" 


they obtain a cycle of length 2“. We see immediately in this case that the cycle 
cannot be of length 2“. As they demonstrate in the same paper, for an z in this 
form (i.e., = 5 (mod 8)) the periods of successive digits from the least significant 
digit are as follows: 1, 1, 2, 4, 8, 16, --- , the least significant digit being always 1, 
the next always 0, the next alternately 0 or 1 the next 0, 0, 1, 1 and so on. 

Now, if we choose p to be of the form 1 (mod 5), say p = 256 = 2°, then the 
effect of multiplying the successive powers of x by 2° and reducing the product 
modulo 2“ is to shift each digit to a position 8 places more significant and to fill 
in the 8 least significant places with zeros. Because of the strict increase by factors 
of 2 in the period of the digits as they become more significant, this must. result 
in a reduction of the length of the maximum cycle by the factor 2°. (No formal 
proof of this maximum appears to have been given previously. ) 

This example should suffice to illustrate the importance of obtaining a set of 
necessary and sufficient conditions on p and x for the maximum cycie-length to be 
achieved, and also of giving a description of the behavior of the system when these 
conditions are not satisfied. 

The only attempts to define formally the restrictions necessary on p and zx 
are those of Leslie and Gower [4] and Certaine [5]. 

Leslie and Gower state that a maximum period of 2*’ distinct random elements 
is obtained subject to: 

(1) Choosing p and 2’ at random; 

(2) Replacing x’ by x the closest number to zx’ such that x = 5 (mod 8); 

(3) Forming successively the numbers px” (mod 2"), r = 1, 2,---. 

The random elements are then all 2“? numbers (mod 2") whose least significant 
binary digits are 10. This result is attributed to a theorem by Euler. 

On examination we find that condition (2) on z is more restrictive than is 
necessary; and, as we have already seen, it is essential that some restriction be 
placed on p. 
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A discussion of the form of z necessary to ensure a cycle of maximum length 
and of the length of this maximum cycle has been given by Certaine [5] for general 
M, but his conclusions on the form of z for the particular case 


M=2 


would again appear to be unnecessarily restrictive. 

Furthermore, no attempt has been made to describe in what way the system is 
affected by a choice of x and p which do not satisfy the conditions required for the 
maximum cycle. 

The system of numbers obtained from equation (1) for M = 2" is fully de- 
scribed in the following section, the proofs of the results being given in Section 3. 


2. Formal Description of Systsm of Numbers Generated by z,., = x-z, (mod 2"); 
p = xo. Under favorable conditions on z and p we obtain a maximum cycle of 2°” 
elements, all of which are distinct. The conditions on z and p to achieve this maxi- 
mum cycle are: 


I xz = +5(mod8), ie. 
xz =5(mod8) or z = 3 (mod 8) 
II: p = 1(mod2), ie. 
p must be odd. 


I and II are necessary and sufficient for the maximum cycle to be obtained. 
Relaxation of these conditions affects the length of the cycle, and in some cases 
causes the process to degenerate to zero. 


2.1. Relaxation of Condition II. If p is even and condition I is satisfied, the 
maximum cycle of 2°” distinct iterates is reduced in length by a factor 2’ where 
this is the highest power of two by which p is divisible, i.e., if p = 2’ (mod 2”**) 
the cycle is of length 2°-**. 

We have already seen an illustration of this effect in the discussion of the system 
described by Taussky and Todd [1]. 


2.2. Relaxation of Condition I. 
(a) If z is even, the maximum number of distinct iterates is k, generated by 
x = 2 (mod 4). In general, if x = 2’ (mod 2’*’) the number of distinct iterates is 


H (where [z] signifies the least integer greater than, or equal to, z). In all cases 


for x even the process degenerates to zero on the Bo element produced. If p 


is even, say p = 2' mod 2'*’, the process degenerates to zero on the Ea 


element produced. 


(b) If x is odd, the maximum cycle is generated, consisting of 2°” distinct 
iterates, if and only if 


x = 3 or 5 (mod 8). 
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For any other odd values of x the length of the cycle is decreased as follows: 
If x = 7 or 9 (mod 16), length of cycle is 2" *; 
If x 


15 or 17 (mod 32), length of cycle is 2*~, ete. 


We may completely specify ail odd integers in this way, with the general result 
that if z = 2’ + 1 (mod 2’*") the length of the cycle obtained is 2°’; j = 2; all itet- 
ates being distinct. 

We may summarize these results as follows. 

If k = 2, 3 the maximum number of distinct elements is k, generated by x = 2 
(mod 4), p = 1 (mod 2), and the process degenerates to zero. 

If k > 3 the maximum cycle is of 2°” distinct elements and is generated by 
x = 3 or 5 (mod 8), again for p = 1 (mod 2). 

Small values of k are, of course, of little practical interest in the generation of 
pseudo-random elements. 

Relaxation of the conditions I and II simultaneously has the effect of combining 
the results of paragraphs (a) and (b). For example, 


p = 2 (mod 4) 
x = 7 (mod t 

will produce a cycle of 2‘~* distinct elements, or, again, 
p = 2 (mod 4) 
x = 2 (mod 5 


will produce (= ) distinct elements, the process degenerating to zero. 


3. Proof of Results of Section 2. Let us first consider condition I. For general x 
we have two possibilities: 


xeven: ie., x = 2m 
x odd: x =2m+1 
ak 

3.1. x even. M=2 


Let N(x) = [k/I(x)| where I(x) is the index of the greatest power of 2 which 
divides x. Then 
I(a*) = k (since I(x”) = mI (x) # 0). 
Hence x” = 0 (mod 2") and number of distinct iterates cannot be greater 
than N(x), for at this stage the process degenerates to zero. 
Further, by the definition of N(a), the process cannot terminate earlier. 


Also, if 0 < m,n < N(x) > a” = a” (mod 2‘), then 


ei 1 (mod 2‘) 


and hence 


m—n= 0. 
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Therefore, the congruence classes of the N(x) + 1 powers of z: 1, x, 2’, --- , 2", 
are distinct and we must obtain N(x) distinct iterates up to degeneration of the 
process (z° not being obtained). N(x) will be a maximum for N(z) = k, ie., 
I(x) = 1. So x = 2 (mod 4) generates k distinct elements, and in general z = 2’ 


(mod 27**) generates [*}. 
L 


3.2.20dd. M=2;2 =2m+1. 
LEMMA. Z-—I1(Z)25 fon Z25 


Proof. If 2‘ S$ Z < 2 then I(Z) Ss i. Thus Z — I(Z) = 2‘ — i which is mono- 
tone increasing in 7 and so 25 fori = 3, ie., for Z = 8. Also if Z = 5, 6,7 then 
I(Z) = 0, 1, 0 respectively. Hence Z — 1(Z) 2 5, all Z = 5 as required. 

Now if z = 2m + 1, the 2°" congruence classes of z in the range (0, 2") form a 
multiplicative group, the order n(x) of the congruence class of z being the least 
integer such that 


z"® =] (mod 2, 


hence, n(x) divides 2"~ {6}. 
Furthermore, the congruence classes of 2, x”, --- , 2” include all powers of 
x and are distinct, for if 


= x” (mod 2") m>n, then 
z”-" = 1 (mod 2‘), 


and m — n must be greater than n(z). 

So we see that the cycles generated by odd integers have length of a power of 
2, the maximum cycle being given by the odd integer of greatest order in the group, 
and we must find the conditions necessary on z for it to have this greatest order. 
That the process must cycle is cbvious, for if 


x’ = 0 (mod 2"), then 
xz =0(mod2*) since x 2 


i. e., x is even, which is not true. 
Now, n(x) is determined by 


(1+ 2m)*' # 1 (mod 2*) 
but (1 + 2m)""" = | (mod 2"), where 1 = I(n(x)) — 
Consider (1 + 2m)" (mod 2'*?) for suitable p. 


If p = 3, then 
1{ (omy (? )} 
rj) 
Thus 


4 l 
(2) (1 + 2m)" = {e (2m)’ (7') (mod 2'**), 
0 


—_ 


IV 


r—I(r)+l 


>I1+5 for r 


IV 
oe | 
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and setting ¢ = (1 + 2m)* — 1, we have 
(+3 
¢ = [2m of. g***(9' a 1)m? + ~ (2! ce 1)(2"* ae 1)m® 


+ a. (2' — 1)(2°* —1)(2' — 3)m'| (mod 2'**), 
Hence, if 1 => 2 
¢ = 2'"im — m? — 2m) (mod 2'**) 
= 2'"'[m(m + 1) — 2m*(m + 1)] (mod 2'**) 
= 2'**m(m + 1) (mod 2'**) 


(for if m odd or even 2m?(m? + 1) = 0 (mod 4)). 
We must now distinguish between 


(i) m 


Il 


0, —1 (mod 4) 
1, 2 (mod 4). 


(ii) m 


(i) (1 + 2m)” — 1 = 0 (mod 2'**) hence, if 1 + 3 = k, x™’* = 1 (mod M)) 
so that n(x) divides M/8. 

(ii) (1 + 2m)" — 1 = O (mod 2’) (for m(m + 1) = 2 (mod 4) 
but (1 + 2m)” — 1 # 0 (mod 2"). Hence, taking k = 1 + 2;k = 1+ 3, we have 
a™'* = 1, 2™"* % 1, (mod M), so that the maximum value of n(x) is M/4, and it 
assumes this value for all x = 1 + 2m, where m = 1, 2 (mod 4), i-e., 


x =3(mod 8), or 
x = 5 (mod 8) 


and only these values. 
Taking p = 4, we again find that 


= m(m + 1)2'*" (mod 2'**) (l = 3) 


and we can immediately determine for what values of x we obtain the next largest 
cycle. 
For this we want 


m(m + 1) = 0 (mod 4), and 
m(m + 1) #0 (mod 8), i.e. 


m 


3,4 (mod 8), giving 
7,9 (mod 16). 


zx 


Similarly, we may extend this result by suitable choice of p to show that in 
general we obtain a maximum cycle of 2°’ distinct elements for 


az = 2’ + 1 (mod 2***). 
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This completely specifies the system for all odd z, for any odd value of z can be 
expressed in the form 


x = 2’ + 1 (mod 2”), j= 2,3, --- 
(x = 1 is, of course, trivial). 

We now consider the effect of different values of p on these results. Any value of 
x which is of the form z = 2’ + 1 (mod 2’); 7 = 2, i.e., any odd value of z, pro- 
duces a cycle of 2*~’ distinct elements having the following characteristics: 

(a) The least significant digit is 1; 

(b) The least significant 7 + 2 digits are of total order 4; 

(c) The order of the (7 + 3)th digit is 2, the order of the (7 + 4)th digit is 4, 
etc.; that is, the order of the successive digits beyond the (j + 2)th have orders 
which increase by the factor 2 as they become increasingly more significant. 

Therefore, it is apparent by inspection that the effect of multiplication of the 
elements by any odd value of p will be to leave these characteristics invariant and 
to unalter the length of the cycle of elements obtained. 

Furthermore, if p is even, say p = 2” (mod 2"**), the r least significant digits 
must become zero and the above characteristics will then be true for the remaining 
k — r digits, i.e., multiplication by an even-valued p has the effect of shifting the 
digits to a position r places more significant, performing a permutation of the digits 
remaining after such a shift, which does not affect their order, and substituting 
zeros for the least significant r digits. This must result in a reduction in the length 
of the cycle of elements by the factor 2”. 

These remarks are easily verified if we consider the effect of multiplication by 
p of the equation (2) for suitable choice of p for the form of x to be studied. 

If, however, x is even, say z = 2’ (mod 2’*"), the successive powers of z have j, 
2j, --- zeros as their least significant j, 27, --- digits, and multiplication by 
any odd-valued p cannot affect the number of distinct iterates before degeneration, 
for it must ensure that in the successive iterates the (j + 1)th, (27 + 1)th--- 
digits are non-zero. 

When »p is even, say p = 2” (mod 2”**), multiplication by p will introduce r 
further zeros in piace of the r least significant non-zero digits and will therefore re- 


duce the number of iterates to E rt. 
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A Note on the Convergence of Alternating 
Direction Methods 


By Milton Lees 


1. Introduction. Convergence of the single-parameter alternating direction 
methods of Douglas, Peaceman, and Rachford [1], [2], [3] has been proved for a wide 
class of elliptic difference equations; see, for example, Birkhoff and Varga [4]. The 
proof consists in showing that a certain matrix, similar to the defining matrix, has 
spectral radius less than one. Since the defining matrices for these methods are 
symmetric only when they are induced by a proper discretization of Laplace’s equa- 
tion in a rectangular region, an estimate for their spectral radii does not imply a 
corresponding norm estimate for their rate of convergence. 

The purpose of this note is to present another proof of the convergence of the 
two basic alternating direction methods which, at the same time, provides a norm 
estimate for their rate of convergence. First, for simplicity, we shall consider 
Laplace’s equation in an arbitrary, bounded lattice region. We shall prove that it is 
possible to select the acceleration parameter so that a suitable norm of the error is 
reduced after each iteration by the amount 1 — ch + 0(h’), where h is the uniform 
node spacing for the lattice, and c > 0 depends only on the minimal eigenvalue 
for the Laplace operator. We shall indicate only briefly the extension of these re- 
sults to elliptic equations with variable coefficients. 

Our proof of the convergence of the alternating direction methods was motivated 
by two considerations; first, the close resemblance of the iteration equations to 
parabolic difference equations, and second, an integral estimate technique for 
proving that, under suitable conditions, the solutions of parabolic differential 
equations decay exponentially. 


2. Formulation of the Problem. Let £ denote the (uniform) lattice of side h 
determined by the nodes (ah, 8h), a and 6 being integers, positive, negative, or 
zero. Denote by §S any finite subset of £. In the usual way [5], we decompose § 
into two disjoint subsets: S, the interior nodes of S, and 4S, the boundary nodes 
of 8S. 

If Q is any subset of £, we denote by C(Q) the linear space of all real-valued 
functions on £ whose support is contained in Q. (In general, the support of a function 
f is the closure of the set {p | f(p) # 0}.) Note that @(S) is finite-dimensional with 
dimension equal to the number of nodes in S. The advantage of thus trivially ex- 
tending functions on Q to all of £ will become apparent later. 

We denote by A the usual “‘five-point”’ Laplace difference operator, i.e., 


Au = Use + Uys, u € C(Q), 


where the subscripts z and y denote forward difference quotients, and ¢ and @ de- 
note backward difference quotients [6]. 


Received July 20, 1961. The work presented in this paper was supported by the AEC Com- 


puting and Applied Mathematics Center, Institute of Mathematical Sciences, New York 
University. 


70 

















CONVERGENCE OF ALTERNATING DIRECTION METHODS 71 


We shall be concerned, first, with the Dirichlet problem for the elliptic difference 


equation 
(1.1) Au=f in 8, 
where f « €(S). More precisely, we seek a function u « @( 8) that satisfies (1.1) 
and assumes arbitrary prescribed values on 98. 

It is easy to see that this problem is equivalent to the problem of solving a linear 
system of algebraic equations of order equal to the number of nodes in S. By means 
of the familiar maximum principle for the difference operator A, one can show that 


this system of linear equations always has a unique solution. In this paper we shall 


be concerned with the two basic, single-parameter alternating direction methods 
for the inversion of this system of equations. 


3. The Alternating Direction Methods. Let w(0) « €(8) be a function that agrees 
with u on 0S, and otherwise arbitrary in S; w(0) will serve as the initial approxi- 
mation to u. Now, for n = 1, we determine two sequences of functions {w(n)} and 
{w*(n)} in C( 8) such that 


(2.1) w(n) = w*(n) = u on 8, 
and 


(a) p [w*(n) — w(n — 1)] = wa(n) + wy(n — 1) —f 
(2.2) 


(b) p [w(n) — w*(n)] = wyg(n) — wyg(n — 1) 


in S. Here p > 0, the so-called acceleration parameter, is to be selected later. These 


equations define what is usually referred to as the Douglas-Rachford alternating 
direction method for solving Au = f. 


In [3] it is shown that (2.1) and (2.2) determine uniquely (alternately) the 
functions w*(n) and w(n) in terms of w(n — 1) and the values of u on 9S. This 
operation taking w(m — 1) into w(n) involves only the inversion of approximately 
2./p tridiagonal matrices, where p is the number of nodes in S. 

Following Douglas and Rachford, we eliminate from (2.2) the auxiliary fune- 
tion w*(n). Adding (2.2a) to (2.2b), we obtain the equation 


(2.3) Dw(n) = wi(n) + wyg(n) — f, 
where we have introduced the notation 
Dw(n) = p '[w(n) — w(n — 1). 

Solving (2.2b) for w*, we get 

w*(n) = w(n) — p Dwy,(n). 
Hence, 

wrs(n) = wWrs(n) — p ADw(n), 
where A is defined as follows: 


Au = Uziyg - 
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Putting the value just obtained for w2,(n) into (2.3), we find that w is a solution 
of the fifth-order difference equation 

(2.4) Dw = Aw — pADw —f, in S, for n=l. 


The other alternating direction method, the Peaceman-Rachford method, is 
very similar to this; it is characterized by the equations 


p {w*(n) — w(n — 1)] = wa(n) + wy(n — 1) —f 
and 
p {w(n) — w*(n)] = wie(n) + wy(n) — f. 


After eliminating from these equations the auxiliary function w*(n), we find that 
w is a solution of the fifth-order difference equation 


(2.5) Dw = 2Aw — pADw — p' ADw — 2f in S, for n2= 1. 


. Our problem, now, is to show that w(n) — u, asn — «. 


4. Preparation. We define on @(S) an inner product (u, v) as follows: 
(u,v) = WD) u(P)o(P). 
PeL 


The norm (v, v)"” induced by this inner product will be denoted by || v ||. We re- 
quire two additional norms for @(S): 


ll» Ia = (I ve I? + |] 96 |)" 
and 


Il > lle = (ll ol]? +o | vag |I?)*”. 
A subscript S will be placed on these quantities whenever we want to indicate that 
the sum is to be extended only over S. 


In [7] it is proved that the minimal eigenvalue \ of the Laplace difference opera- 
tor A, relative to S, can be characterized as follows: 


ll» [lis 


A= UAE 
Oxvee(s) || v ||? 





Since || » || 1,s < || v || 1, it follows immediately that 
(3.1) Allo |? S lolly, 


for all v e C(S). 
On the strength of this inequality, we now prove the 
Lemma 1. Any function v « C(S) satisfies the inequality 


(3.2) ll v |x” > Aoll v |] 2, 
where 

r 
ani NO TF On 


Proof. We have from (3.1) and the definition of || v || » that 


Il» |’ — Allo ls’ = —Apl| ves ||”. 
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To this inequality we apply the elementary inequality [8] 
ll v0 ||? S 2\) oI x, 
valid for all v « €(S), and we find that 
I] o llr — Allo jr > —2Ap'h* |] o |] Y’, 

which is equivalent to (3.2). 

5. Convergence of the Alternating Direction Methods. Let u « @(8) be the 
solution of the Dirichlet problem for Au = f in S, and let w(n) be the correspond- 
ing solution of the Douglas-Rachford equation (2.4). Then, by linearity, the error 


function v(n) = u — w(n) belongs to the linear space €(S), and satisfies, form = 1, 
the difference equation 


(4.1) Pv = Dv — Av + p ADv = 0 


in S. Let o(n) = (1 + pu) ”, where uw > O. One verifies immediately that ¢ is a 
solution of the difference equation 


(4.2) Do = —ud. 


THEOREM 1. Jf » S Xo, then the error function v(n) for the Douglas-Rachford 
method satisfies the inequalities 


| v(m) |] « S o(m)|] v(0) | « (i = 1,2) 
Proof. We set 


(4.3) M = 2p do *(&) (v(), Pv(é)). 


Since v(£) « C(S), we see that M = 0. 
Introduce the function z(¢) « €(.S) by means of the formula 


(4.4) v(&) = o()z(é). 
Then we have that 
Dv = @Dz + zDo — pD¢o-Dz 
which, in view of (4.2) becomes 
Dv = ¢laDz — uz], 


where a = ¢ (1) > 0. 
From this and (4.1) it follows that 


(4.5) Pv = ¢faDz — pz — Az + ap ADz — p'yAz] = olaDz + Bz], 


where B is defined in the obvious way. 
Because of (4.5), (4.3) becomes 


(4.8) a »d2 (e(£), aDe(t) + Be(8)). 





74 MILTON LEES 


VY e now recall from [8] the following quadratic difference identities 


(4.7) 2 (z, Dz) = D\| zl’ + ell Dz |I, 
(4.8) 2 (2, ADz) = D|| 249 ||’ + ol] Dees ||’, 
(4.9) (2, Az) = —|[z| 1, 

and 

(4.10) (2, Az) = || 22g I’. 


Using these identities in (4.6), we obtain, after some slight rearrangement, the 
identity 


(4.11) M= pd {aD|| z ll’ + 2(|| 2 lly’ — w || 2 lle) + || Dz |x’ + wal] Dz |\,'}. 
Since || z ||," = oll z ||.’, by Lemma 1, and u < Xo, we have from (4.11) that 


M2p > aD|| z ||2' = al] z(”) |l2"— al] 2(0) |-’. 


But, since M = 0, we obtain from this and (4.4) that 
|| v(m) |lz = o(n)|| 2(m) lle S o(n)|| 2(0) 2 = O(n) |] v(0) lle, 


which proves the second inequality of the theorem. 
Next, we set 


(4.12) N=p > ¢ *(€)|] Po(é) |Is°. 
Then, by (4.4) and (4.5), (4.12) becomes 


N= pd || aDz(€) + Be(é) ||s°. 
But 
|| aDz + Bz ||s’= a’ || Dz ||s'+ || Bz ||s’"+ 2a(Dz, Bz)s = a’|| Dz ||* + 2a(Dz, Bz), 
where we have dropped the subscript S, because Dz « C(S). Hence, 


(4.13) N = pd {a"|| Dz ||’ + 2a(Dz, Bz). 
t=1 
Using the identities (4.7)—(4.10) again, we obtain from (4.13) the inequality 
N2p 2, {a°|| Dz ||’ + aD\| z ||’ — waD|| z |\2"+ ap|| De |\:}. 


Therefore, by adding a combination of the above inequalities, 


N + uM = pa > Di z llr. 
Again, N + uM = 0, so that 
| 2(m) |], S |] 2(0) |, 
which, in view of (4.4), implies that 
|| v(m) |], S o(m)|| v(0) Il. 
This completes the proof of the theorem. 
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Since pp > 0, Theorem 1 implies that the Douglas-Rachford method is con- 
vergent, but we also have the stronger result: 

THEoreM 2. If u = \oand p = h/+/2X, then the error function v(n) for the Douglas- 
Rachford alternating direction method satisfies the inequalities 


(4.14) | o(n) ||; s (1 + = n) 1 »(0) I, (¢ = 1,2). 


Proof. The inequality (4.14) follows from Theorem 1, by minimizing ¢(n) as 
a function of p. 


In an analogous fashion, we can prove 
THeorem 3. If u = 2do(1 — dop) * and p = h/+/2X, then the error function v(n) 
for the Peaceman-Rachford alternating direction method satisfies the inequality 


o t= vB ion 
|e(n) fh = ($= Vey" | (0) Is 


Note that the Peaceman-Rachford method has a smaller convergence factor. 
The preceding results can be extended to more general elliptic difference equa- 
tions, e.g., to those that come from elliptic differential equations of the form 


a s du\ 


in a bounded, open set @. The functions a and b must belong to c'(@). The main 
point that we wish to emphasize is that, because the inner product (wu, v) is defined 
as a sum over all of £, one must be able to extend a and b in a continuously differ- 
entiable way to a slightly larger set w > ©. This can be carried out, for instance, by 
means of the Whitney extension theorem [9], in a form given by Hérmander [10]. 

One can also show that similar results can be proved for the Douglas-Rachford 
method as applied to elliptic equations in n independent variables, although this 
is not the case with the Peaceman-Rachford method [8]. 
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Calculation of 7 to 100,000 Decimals 
By Daniel Shanks and John W. Wrench, Jr. 


1. Introduction. The following comparison of the previous calculations of x per- 
formed on electronic computers shows the rapid increase in computational speeds 
which has taken place. 

















Author Machine Date | Precision Time 
Reitwiesner {1} ENIAC 1949 | 2037D 70 hours 
Nicholson & Jeenel [2] NORC | 1954| 3089D | 13 min. 
Felton {3} Pegasus 1958 | 10000D 33 hours 
Genuys [4] IBM 704 | 1958 | 10000D | 100 min. 
Unpublished [5] IBM 704 1959 | 16167D 4.3 hours 





All these computations, except Felton’s, used Machin’s formula: 
(1) x = 16 tan} — 4 tan’ ss. 


Other things being equal, that is, assuming the use of the same machine and the 
same program, an increase in precision by a factor f requires f times as much 
memory, and f° times as much machine time. For example, a hypothetical com- 
putation of x to 100,000D using Genuys’ program would require 167 hours on an 
IBM 704 system and more than 38,000 words of core memory. However, since the 
latter is not available, the program would require modification, and this would ex- 
tend the running time. Further, since the probability of a machine error would be 
more than 100 times that during Genuys’ computation, prudence would require 
still other program modifications, and, therefore, still more machine time. 


2. A New Program. We discuss here a computation of x to more than 100,000D, 
which required 8 hours 43 minutes on an IBM 7090 system. This increase in speed 
by a factor of about 20 is largely due to the increased speed of the 7090 (it is about 
7 times as fast as a 704), but substantial gains were also obtained by programming 
changes. 

The formula we used, namely, 


(2) we = 24tan™ +8 tan’ sr +4 tan 355; 


is due to Stérmer [6], and has not been previously used in high-precision computa- 
tion. The computation time breaks down as follows: 
8 tan #y: 3 hours 7 minutes 


1 


4tan $y: 2 hours 20 minutes 


24 tan’ 4: 2 hours 34 minutes 
Received September 7, 1961. 
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Conversion (binary-to-decimal ): 42 minutes. 


To obtain these favorable times, two devices were used. 
a) Instead of evaluating the Taylor series 





al = & (-1)‘Am 
(3) A tan er & (2k + 1)m2@+ 


term by term, one may compute two terms at a time by using 


“s =. Am|(4k + 3)m® — (4k + 1)] 
(4) Ata oS" 2, (16H + 16k + 3)m**) 
This substitutes 2 (multi-precision) divisions, 1 multiplication, and 1 addition for 
4 divisions, 1 subtraction, and 1 addition, and this eliminates 27% of the computa- 
tion time. Further improvement in this direction, i.e., three terms at a time, is not 
possible here, since the divisors in (4), m* and (16k° + 16k + 3), already tend to 
fill a 7090 word, whose maximum numerical size is 2” — 1. 

b) Since the 7090 is a binary machine, the (multi-precision) division by m‘ and 
the multiplication by [(4k + 3)m? — (4k + 1)] may be replaced by a simple shift- 
ing operation for the case m = 8 in 24 tan™’ 4. Had this not been done the 2 hours 
34 minutes listed above would be, instead, 6 hours 7 minutes with the double-term 
formula, (4), or 8 hours 21 minutes with (3). 

The program using (4) requires three blocks of storage, one for the current 
Am/m*“*» | one for the current term, and one for the partial sum. Since a 7090 word 
is equivalent to 10.536 decimal digits, the working storage for an N-place 7 is 
3N/10.536. Therefore, a core memory of 32,768 easily suffices for a computation 
to more than 100,000D. 





3. The Check. Inasmuch as such a computation requires billions of (arithmetical) 
operations, it is clear that a check is necessary. This was obtained with Gauss’s 
formula [7]: 


(5) x = 48 tan’ + + 32 tan” 2 — 20 tan’ shy. 
During the computation, using (2), the numbers 


A = 8tan 
and 
B = 8 tan’ 3 + 4 tan‘ sts 


were written on tape. At the start of the check A and B were read into memory 
and 9A — 5B = 32 tan” x — 20 tan $5 was computed (in 1 second). To this 
difference was added the new number 48 tan” 7s, which was computed by (4) in 
4 hours 22 minutes. The check, therefore, takes less time than the original run 
(8 hours 1 minute), and is perfectly valid, since an error in any of the four are- 
tangents will lead to a discrepancy between the results of (2) and (5). 


4. The Result. The run and the check were both made on July 29, 1961, and 
such a discrepancy did in fact occur. The two values of z, in binary form, were com- 
pared by the machine in } second, and were found to agree to only 234,848 bits. This 
is equivalent to 70,695 decimal places. Subsequently the error was isolated. It 
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was found to have occurred in the computation of 24 tan’ 4. The second value of 
=, computed using (5), was therefore correct throughout. When 24 tan } was 
recomputed, the two values of x agreed up to the last word. This comprised 333,075 
bits, or 100,265 decimal places. The first 100,000 of these are given here. 

The computation was of such a character that it was known a priori that the 
term involving tan 3}5 contains most of the round-off and truncation error, and 
consequently that we have the inequalities: 


x computed by (2) < x < «x computed by (5). 


Thus the check and the computation give upper and lower bounds respectively, 
and, to the extent that they agree, x is determined absolutely. 

Care has been taken in the output routines—in writing and printing—and, since 
the reproduction here is photographic, we believe that this value is entirely free 
from error. 


5. A Million Decimals? Can x be computed to 1,000,000 decimals with the com- 
puters of today? From the remarks in the first section we see that the program which 
we have described would require times of the order of months. But since the memory 
of a 7090 is too small, by a factor of ten, a modified program, which writes and reads 
partial results, would take longer still. One would really want a computer 100 times 
as fast, 100 times as reliable, and with a memory 10 times as large. No such machine 
now exists. 

There are, of course, many other formulas similar to (1), (2), and (5), and other 
programming devices are also possible, but it seems unlikely that any such modifica- 
tion can lead to more than a rather small improvement. 

Are there entirely different procedures? This is, of course, possible. We cite the 
following: compute 1/x and then take its reciprocal. This sounds fantastic, but, 
in fact, it can be faster than the use of equation (2). One can compute 1/z by 
Ramanujan’s formula [8}: 


Sn Te 
w 4\882 882? 2 4 882° 2-4 4°. 8? ; 

The first factors here are given by (—1)* (1123 + 21460k). A binary value of 1/z 
equivalent to 100,000D, can be computed on a 7090 using equation (6) in 6 hours 
instead of the 8 hours required for the application of equation (2).* To reciprocate 
this value of 1/x would take about 1 hour. Thus, we can reduce the time required 
by (2) by an hour. But unfortunately we lose our overlapping check, and, in any 
case, this small gain is quite inadequate for the present question. 

One could hope for a theoretical approach to this question of optimization—a 
theory of the “depth” of numbers—but no such theory now exists. One can guess 
that e is not as “‘deep”’ as x,f but try to prove it! 

Such a theory would, of course, take years to develop. In the meantime—say, 
in 5 to 7 years—such a computer as we suggested above (100 times as fast, 100 
times as reliable,and with 10 times the memory) will, no doubt, become a reality. 
At that time a computation of x to 1,000,000D will not be difficult. 





* We have computed 1/x by (6) to over 5000D in less than a minute. 
t We have computed e on a 7090 to 100,265D by the obvious program. This takes 2.5 
hours instead of the 8-hour run for + by (2). 
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TECHNICAL NOTES AND SHORT PAPERS 
Characteristic Exponents of Mathieu Functions 


By T. Tamir 


The study of the properties and solutions of the Mathieu second-order differ- 
ential equation has been very extensive in the last decades, due to the special 
interest presented by physical problems involving periodic media and problems 
separable in elliptic coordinate systems, [1], [2]. The purely periodic solutions of 
the Mathieu equation have been extensively tabulated with high accuracies over 
large ranges of the parameters involved [3], [4]. 

On the other hand, numerical tables for the non-periodic Floquet-type solutions 
are scarce [5], [6], [7]; in addition, the increments between two successive tabulated 
values are relatively large for some of these tables and the accuracy is relatively 
low in the other tables. Consequently, the applicability of the available data is 
rather limited. 

The regions of “stable” solutions of the Mathieu equation are of particular 
interest in applications involving the wave equation with periodic boundary con- 
ditions, since these regions correspond to the propagating waves in the correspond- 
ing media. The first three regions of stable solutions were tabulated with this type 
of problem in mind, the numerical results being given in Tables I, I, and ITI. 

The canonical form of the Mathieu equation considered is given by 


(1) y” + (p — 2q cos 2z)y = 0, 


with p and qg real. Every non-periodic solution of (1) is a linear combination of two 
linearly independent Floquet-type solutions of the form 


(2) y = e”"P(z), 


where P(z) is a periodic function of z, and v is the so-called characteristic exponent. 
This exponent is a function of the parameters p and gq only. When » is real, the 
solutions (2) are called “stable” since they are uniformly bounded for all real z; 
then v = f(p,q) may be chosen to be single-valued. (Note: In the tables and curves 
shown here, v is not reduced to 0 < v S 1 as is usually done, but its actual value 
is taken, thus defining the particular stable region it represents. ) 

It is shown in the references that all values of v are solutions of the continued 
fraction equation 





2S ee 1| 1 | 1 | 
3 D ed —in, —ele “ee — — = — = — 
(3) . D, | De | Ds 2 Ds | D_» | D_s 
where 

2 
(4) = OO te . m<tate3---). 
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For computation purposes, it is convenient to define 














(5) z=v—p 

so that equation (3) will take the form 
ORE TCS. td RE Fae eae ee 
4i+y+2 |82+ 7) +2 |4n(n + v) +2 

(6) Wi oe q a ari ie ¢ F 
4i-vy+2 |82-—y7+2 |4n(n — vy) +z 


(n = 1, 2, 3, +++) 


The continued fractions are convergent since they satisfy the convergence test 
[2] which requires that | D, | = 2 for n > N, where N is a finite integer. 

To find the roots of equation (6), it is noted that this equation is already ex- 
pressed in a suitable form for solution by means of an iteration process, with z 
regarded as the unknown. Taking g and » as variable parameters, such an iteration 
procedure was programmed on a computer which obtained p with an accuracy of 
10°. As is proved in the Appendix, the iteration process converges in an alternating 
fashion, i.e., the exact result always lies between the values obtained by two suc- 
cessive iterations; hence, the iteration process must be carried out only until the 
difference between two successive results is less than the maximum desired devia- 
tion. 

It was also noted that the number of iterations required to solve equation (6) 
was extremely large for a large range of parameters if a simple iteration procedure 
was used. This was due to a very slow convergence of the iteration process in some 
regions of g and v. To improve this situation, the computer was programmed to 
differentiate between two cases: 

a) “Fast” convergence: x,4; contained within the interval (z, , 2.42), 

b) “Slow” convergence: x4; outside the interval (7, , 2.42), where 2,4; is 
the result of the nth iteration, using z, as the trial value; z.,, is the arithmetic 
mean of 2,4; and x, . 

After using the test to determine whether the convergence is “fast’’ or “slow,” 
as defined above, the computer used either 2,42 or 2,42, respectively, as the next 
trial value for the n + 2 iteration. This procedure reduced the number of itera- 
tions required from more than a thousand to less than twenty for the most slowly 
converging cases. 

The results thus obtained are given in the appended tables, while Figure 1 
ese these values plotted in the p — q plane. Note that, for greater clarity, 

= + ~/'p' rather than p is plotted in the figure. Thus p = P ifP20;p=—-P 
: P <0. 


Appendix. It is shown below that an iteration process based on equation (6), 
if converging, will do so in an alternating manner; i.e., the value of the exact solu- 
tion always lies between the values of the results of any two successive iterations. 

The iteration process may be written in the form 


= Ro* (an) + Ro (2n) 
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FIRST REGION SECOND REGION — THIRO REGION 
af q ice Se tease, 






































aon. 
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v=3.0 
20 v=20 d 
| v=2.5 
v=0.5 | 
| 
| 
. LLL 
=10 a 1) 5 10 5 1) 2c 2.5 30 


Fic. 1—Curves of q vs P for constant values of ». The values of » pertaining to each 
curve are given directly by the interceptions with the P axis. p = P*?if P => 0;p = — P* if P <0. 


where 
eer rps Ah matecm@ied Mf swe oye! Rez | 
Ro (an) 7 Qo* In 4 bot lay*2, + b= (an, +b ae 
Then 
+ u 1 Pit e315 1 | yr 
Rate) = Ge Ft [abate + Ee 


is a “remainder” function or the “tail’’ of a continued fraction after the first m 
terms have been omitted. 
Disregarding the + sign, we have 


1 
Go Ln + bo — Ri(an)” 
Differentiating with respect to x, , we obtain 
ad a — Ry’ (xp) 
[a1 %, + bi — Ri(x,)]? 
Noting that R,,(z,) has the same functional form as Ro(z,), one finds 


Ro’ (xn) = —laRe + a;(RoR;)’ + a2(RoRiR2)° a ee |. 


Ro(2n) = 





Ro’ (an) = 





= —R¢ [ap ~— Ry'(2,)]. 


. 1 : 
In the case considered, am = e or 1; hence am > 0. Rnu(t,) ~ 0 as m— @ 
so that Ro'(z,) has a finite negative value for any finite x, . Consequently, the 
derivative of z,4: is finite and negative so that the iteration process based on 


equation (6), if convergent, will converge in an alternating fashion. 
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CHARACTERISTIC EXPONENTS OF MATHIEU FUNCTIONS 








TZg81'6 
LILLI'6 
SE89T 6 
TS6ST *6 
SOOST '6 
S90F1'6 
60TET 6 
LPIZI'6 
Z8IIl 6 
SIZOI 6 
09260 '6 
9T€80'6 
688206 
S8P90 6 
Z1990'°6 
PLLEO 6 
626806 
VEZEO 6 
9F9Z0'6 
£2610'°6 
ZLETO 6 
&0600 *6 
TZS00 6 
88200 6 
T9000 °6 


LU 





SPESh 8 
PEEFE 8 
6EPEF'S 
Z69ZF 8 
Z80ZF 8 
809TF'S 
OLZIPF'S 
L90IF'8 


6% 





L662 '°8 
16192 °8 
196228 
062618 
Z2291'°8 
9SZEI'8 
£06018 
€Z180'8 
&2990'8 
LO1€0'8 
$8200°8 
20986" 2 
82996" 2 
00276 ° 2 
92626 ° 2 
POPT6' 2 
98668" 2 
02288 2 
L0928° 2 
9F998 "2 
PESSS LZ 
LISS L 
699F8 "2 
€6EF8 2 
LOPS" 2 


8°S 





9866L' 2 
O9TOL 2 
SEPel L 
SE889  Z 
Teego "2 
6L619 2 
9928S" 2 
68999" 2 
8LL29°L 
62009 ° 2 
OFPLE'L 
TEOSP 2 
L8LEPL 
STLOP 2 
9T88E 2 
680282 
vEesse LZ 
SPIEL 
SE6ZE 
Z88TE" 
26608" 
92Z0€° 
LIL6Z° 
81862" 
08062" 


LG 


f+ Py Py ly Py b+ I~ 





PIES LZ 
L818 2 
GSShs 2 
9gF0% "2 
9OS9T* LZ 
GIZZ1'L 
28060° 2 
LE990° 2 
69€Z0 2 
88266 °9 
00496 °9 
GOLE6 9 
L0Z16°9 
90688 °9 
00898 '°9 
06878 °9 
€LTEs8'9 
9F9T8'°9 
L008 °9 
g¢l6l°9 
G8T82°9 
G6ELL°9 
€8292°9 
8PE9L'9 
Z8092°9 


9°% 








9F968 "9 
sors 9 
€9262°9 
96092 °9 
86S0L°9 
£8299 9 
19129 °9 
TPZ89 9 
PEShS 9 
crore '9 
SLLLP'9 
SELFr'9 
SZ61h 9 
6E86E "9 
O869E "9 
SPSrE 9 
T86ZE 9 
SESE °9 
64162 9 
PLP8Z ‘9 
€0FLZ°9 
€ES9S 9 
09892 9 
Z8ESS 9 
$6092 "9 


c’% 


a 





CO86F 9 
Z8IPP'9 
6F98E "9 
P8ZEE "9 
90T8Z'9 
ZETES 9 
SLE8T'9 
9S8ET "9 
08960 °9 
6g¢s0'9 
66210°9 
L0E86°¢ 
£8096 °¢ 
62126 '°S 
Ehr68'S 
0Z0L8'°¢ 
9S8P8'¢ 
Cr6Z8'°¢S 
O8Z18°S 
Po86L°¢ 
Z998L°S 
G69LL°S 
0S692°¢ 
1Zvol'¢ 
GOT9L'S 





v'S 


OfTST'9 
98980 °9 
£81Z0°9 
LP6S6'°S 
20668 °¢ 
ZL0P8 ¢ 
T8P82°¢ 
osTeL’¢ 
86089 °¢ 
prEtY 'S 
66889 °¢ 
GLLS'¢ 
LL609°¢ 
LOSLY'S 
COthr's 
OFSIh'S 
ZP06E GS 
Pr89E¢ 
CPGPE'S 
ZEEE S 
9261'S 
T6808 
ZS00€ °S 
89P6Z'S 
LI16Z°¢ 














ZZOLS'S 
0962'S 
1802'S 
00849 °¢ 
G2929°9 
ZELOS'¢ 
980Fr'S 
06S28°S 
9hle's 
60992 "¢ 
SE10Z°¢ 
ZHOST'S 
SPE0T *S 
9090'S 
902Z0°S 
£9186 'F 
O&296' F 
$6086" F 
1£806'F 
0688" 
6LEL8'F 
SE198'F 
88I¢8'F 
PESES'F 
OFS 'F 


G'S 








96089 '¢ 
gs96¢'¢ 
ZOETS 
L96Zh°¢ 
OILPE*S 
09992 °¢ 
6FS81°¢ 
SIZOI'¢ 
660£0°¢ 
LPLS6'F 
60288 °F 
LE0Z8 °F 
98Z9L'°F 
L0002 °F 
SPLPO FP 
94009 'F 
1Z6SS'P 
8LEZS FP 
OOF6P' F 
C69 fF 
OO00SF 
O8PEF FP 
P9ESFE FP 
C6SIh FP 
LYLIP'F 


POETO'S 
£0Sz¢'¢ 
Z89EF 'S 
LESte'¢ 
9F09Z 'S 
LOZLI'S 
TPS80°S 
£6866 ' F 
LPEl6'F 
GE6Z8 ' F 
8L9PL Ff 
61999 F 
£6289 Ff 
6EZ19 F 
ZOOFP  F 
OfI Le F 
£2908 'F 
C89 F 
1ZZ61 °F 
Steril F 
0600T °F 
0€S90 °F 
90280 °F 
89910 F 
9TFO0'F 
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No absolute convergence proof was found for the iteration process itself; on the 
other hand, no convergence instabilities were experienced in the actual computation 
for the range of parameters that was tabulated. 
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An Algorithm for Solving Certain Polynomial 
Equations with Coefficient Parameters 


By Robert D. Stalley 


1. Introduction and General Method. We describe a storage-saving procedure 
that can be used in real time simulation or other situations in which the roots of 
an equation of the form 


lA 
IIA 
IIA 
IIA 


€o ej Cn ; 


(1) »» F(x, t2, wtb: » tm) y =0, m2 2, 0 


e; integral (7 = 0,1,--- ,”), 
must be obtained within severe time limits, i.e., from storage, for values of the 
coefficient parameter m-tuple (2, , %2,-°-- , 2m) from some given set. The expo- 
nents e; are defined so as to be monotonically increasing rathe than strictly in- 
creasing for later notational convenience. 

Suppose there exist relations 


(2) Us = 9;:(%1,%2,°°* > Tm); (¢=put1,---,m), 252435 ™, 
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which give the reduction 

F(a, 22, o2 » Zu) 
(3) = g(X%1, 22, °°* y Lm)fj( Up, Upsr, *** » Um\A (ay, te, --+ , Sm), 


G = 0,1, --- ,n). 
Then, if we substitute from (3) into (1), setting 


(4) yh(a ,%2,°-* , 2m) = BV, 


and divide out g(x: , t2, --- , 2m), we obtain 
(5) Do Silty 5 Ups, “ae » Um )v = @. 
j=0 


Now, if we store the roots of equation (5) for all values of the (m — u» + 1)-tuple 
(Uy , Upt1 » *** ,» Um) Within its range, we can find the roots of equation (1) for a 
given value of the m-tuple (2; , 22, --- , Z») within its domain of substitution by 
computing (ty, , Usi1,°°* , Um) from (2), obtaining from storage the roots of (5), 
and computing the roots of (1) from (4). Obvious separate considerations must 
apply if the value of h(a, , 22, --- , 2m) Or g(%1, Z2,--* , Zm) is zero or if any of 
the above functions are undefined for the given value of (2, , 22, --- , Zm)- 

Basically we have replaced m-parameter storage by (m — uw + 1)-parameter 
storage and a computation which is ordinarily relatively simple. Rather than 
storing solutions point-wise in m-dimensional space, we store them by (u — 1)- 
dimensional manifolds, namely those given by equations (2) taken simultaneously. 
A side calculation, using equation (4), is needed to obtain the roots for a particular 
point of the manifold. 

This reduction is especially interesting when 4 = m, for then we need only single 
parameter storage. Note that if m = 2 we must have » = m. 

We consider two equations for which this reduction can be made. In both equa- 
tions wp = 2. 


2. An Equation of Power-Type. Let equation (1) be 


(6) »» oT] x) y? =0, a; #0, (j= 0,1,---,n), 

where 

(7) a + > ra, =k+re;, (t2,173,°°*,%m,7, k constant; 
@,@2,°** , @» constant for each value of j). 


To reduce equation (6), we first use (7) to obtain 


m m m 
TI] af* = a'ci TT 27" T] 27 = 2 {i (ai"'x" ai". 


t=1 i=2 t=2 t=2 
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Then, after dividing out x;‘, equation (6) becomes 
p oI us')o" = 0, 
j=0 t=2 

where 


—r; 


; = 2, 2, (¢ = 2,3, --- ,m), 


and 
Tr 
v= NY. 


This completes the reduction. 

Homogeneous polynomial equations, with y as one of the variables, are special 
cases where r2 = 73 = --- = rT, = 1, r = —1, and a; is a non-negative integer, 
(¢ = 1,2,---,m). 

Now we look at the case m = 2. Set x, = w, m2 = 2, a; = a, a2 = B, fo = 8, 
r = t, and uw = b. The equation 


(8) 2d, anw*z'y’ = 0, (a + 58 =k + te;), 
3 
is reduced to the equation 
(9) > a'r = 0, 
j=—0 
where 
(10) b=w“z 
and 
(11) v= wy. 


Thus, if the roots of (9) are stored for all of the values of 6 in its range, we have 
the following algorithm for solving equation (8) for given values of w and z: 
I. Compute w’ and w’; 

II. Compute b from (10); 

III. Look up the roots of (9); 

IV. Compute y from (11). 
The algorithm is especially simple if s and ¢ are integers. We have combined storage 
lookup and computation to obtain an algorithm for solving (8) which involves 
little computation and only single parameter storage. The writer once needed a 
real time series solution of an equation of the form of (8) where e, = 6, s = 2, 
t = —4, and where, of course, the remaining arbitrary constants were equal to 
specified values [1]. 

This equation of power type, equation (6) with condition (7), has many general- 
izations which are reducible in the same way. The functions F; may be quotients 
of polynomials whose terms are of the form of the coefficients of equation (6), 
where e; is constant for all terms of F; , where k and r are constant for each poly- 
nomial, and where the values of k and r for any numerator minus the corresponding 

















AN ALGORITHM FOR SOLVING CERTAIN POLYNOMIAL EQUATIONS 109 
values of k and r for the corresponding denominator are constant. Mereover, the 


terms may contain other factors and the parameters may themselves be functions 
of other parameters. 


3. An Equation of Exponential Type. Let equation (1) be 


(12) La [I ai*)y" = Q, a; * 0, (j = 0,1,---,n), 
where 
(13) al] a;’ = ke, (r2,13,°** , Tm, ¢, & constant; 

@ , @,°-* , & constant for each value of j), c ~ 0,k ¥ 0. 


To reduce equation (12), we first use (13) to obtain 
Pat — ref cf] at = en(fL are, 
Then, after dividing out k*', equation (12) becomes 
> aj Iai)" = 0, 
j=l t=—2 


where 


“uz =—%— TN, (¢ = 2,3,---,m), 
and 


71 


This completes the reduction. 


There are specializations and generalizations of this equation of exponential 
type analogous to those for the equation of power type. 


4. Remarks. The restriction that e; be integral, (j = 0,1, --- , m), is not used 
and hence may be removed throughout this paper. With the customary interpreta- 
tion of notation and statements, we may extend this paper to cover the trivial 
case » = m + 1, and consequently the trivial case m = 1. 

An open problem is the enumeration of all examples, from various specified 
classes, for which the method of this paper is applicable. Another open problem, 
perhaps more interesting, is that of solving approximately polynomial equations 
with coefficient parameters by finding and solving approximating polynomial 
equations of a type for which the method of this paper is applicable. The generality 
of equation (6) with condition (7) and equation (12) with condition (13) make them 
promising approximating polynomial equations. 
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Corvallis, Oregon 
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A Remark Concerning the Solution of the 
Dirichlet Problem by Finite Differences 


By Bernard Epstein 


In the classic paper [1] it is proven that the mesh-functions obtained by solving 
a discrete analogue of the Dirichlet problem converge, as the mesh-width approaches 
zero, to a harmonic function which solves the Dirichlet problem in a somewhat 
generalized sense. More precisely, let the function f be continuously differentiable 
in a bounded plane domain G and continuous in the closure G of G, and let the 
Dirichlet integral 


(1) p(s) = ff (2 + 5A ae dy 


be finite. Then the finite-difference method presented in the aforementioned paper 
is shown (under suitable assumptions concerning the smoothness of the boundary ) 
to furnish a function u harmonic in G whose Dirichlet integral D(u) is finite (and 
not greater than D(f)), and it is shown that u agrees with f on the boundary dG 
in the sense that for all sufficiently small values of « the inequality 


(2) If (u —f)’ dxdy < yé 


holds, where S, denotes the portion of G consisting of those points whose distance 
from 0G is less than « and y denotes some positive number independent of «. The 
authors indicate, without supplying the details, that u agrees with f on 0G in the 
more elementary sense that the function u — f, initially defined only in G, becomes 
continuous in G if defined to vanish everywhere on 0G. In this brief note we present 
a proof of this fact, thus showing, without reference to any other method of treat- 
ing the Dirichlet problem, that the finite-difference method provides an existence 
proof for the “conventional” formulation of the Dirichlet problem as well as an 
effective procedure for computing the solution. 

The proof is accomplished by establishing two lemmas. The first of these is a 
strengthened version of the inequality (2). 


(3) Lema 1. lime” If (u — f)* dxdy = 0. 
e>0 8, 
Proof. First we establish (3) in the particular case that G is the unit disc, 
e+y <1. Let v = u —f and let two points with polar coordinates (r, 6), 
(R, 6), r< R <1, be selected. Then 


R R 
(4) oR, 0) — o(+,0) = | (0,0) do = | v4(0,0)0'*-0°** dp. 
Applying the Schwarz inequality and assuming that r > 3, we obtain 


(5) [v(R,6) — v(r,@)|° < ({" Upp i\([ de) < 2(1 —r) vp p dp. 
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- 


Integrating over @ and setting « = 1 — r, we obtain 
2r 

[ [v(R, 0) — v(r,6)}° do < 2e I v, dx dy < 2 I [v2 + v,2] dx dy 
0 8. 8, 


< 2D(v). 
Since D(f) and D(u) exist, D(v) also exists, so the left side of (6) must approach 


(6) 


2r 
zero with ¢. It then follows, by a familiar argument, that lim [ v'(r, 0) d@ 
«+0 0 


exists. If this limit, which we denote by C, were positive, the inequality 
| u(r, 0) dé > 4C 
0 


would hold for all sufficiently small e. Integrating this inequality, we would obtain, 
in contradiction with (2), 


1 
(7) I v'(r, 0) rdr dé > lo rdr = 1 oe —e). 
8, 2 l—e 4 


Having thus established that C = 0, we return to (6), let R — 1 (keeping r mo- 
mentarily fixed), and thus obtain, by an obvious application of the triangle inequal- 
ity, 
(8) | v'(r, 0) d@ < 2D,(v) 

0 


where D,(v) = I [v. + v,7] dx dy. Integrating and noting that D.(v) is a de- 
8. 

creasing function of r(= 1 — e), we obtain 

1 


(9) I v dx dy < 2D,(v) pdp = ¢€(2 — €)D,(v). 
8. l—e 


Since D.(v) approaches zero with «, (9) implies (3). 

For a more general domain with sufficiently smooth boundary it is readily 
seen that the same argument may be applied by introducing a coordinate system 
whose coordinate curves are the normals to the boundary and their orthogonal 
trajectories. 

Before stating the second lemma, we introduce the following terminology 
(cf. [2, p. 481]). A D-function (defined in a domain G) is one that is continuously 
differentiable in G and has compact support, i.e., there exists « > 0 such that the 
function vanishes throughout S, . A D-function, say v, is one which is continuously 
differentiable in G and can be approximated by D-functions in the sense that there 
exists a sequence {v,} of D-functions such that 


(10) lim I (v — v,)* dx dy = 0, lim D(v — »,) = 0. 
n>o G no 
LemMaA 2. o(= u — f) isa D-function. 
Proof. As in the proof of the previous lemma, it suffices to confine attention 
to the case that G is the unit disc. For each positive integer n we define the func- 
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tions g, and h, as follows: 


1 jsceet=- 
n 
(11) gn = jos nx (1 ~r — 4), fut eset 2) A wt ~~. 
n n 2n 
1 
0 : rel-> 





Then the functions v, = vg, are obviously D-functions, and we establish the two 
parts of (10) by the following arguments. 


If (v — v,)* dx dy = If (vh,)” dx dy 


(a) 
= # (vh,)° dx dy Ss I. v dx dy — 0; 
D(v Naa v2) = D(vh,) = Dyn (vhn) & Dyjn(vhn) 
(b) + If [vz hn — URn,z)” + (vy hn — vhay)”] dx dy 


1/n 


er Il [An(vs + vy) + (hie + hin y)] dx dy. 


1/» 


Since h,’ S 1 and hi. + hi, S n'a’, we obtain 
D(v — vn) S 2Dyn(v) + 2n?x? If v dx dy. 
Si s/n 


Dy,(v) approaches zero with increasing n, and the same is true of the remaining 
term, by Lemma 1. Thus the present lemma is proven. 

The desired result concerning the boundary behavior of u now follows imme- 
diately from the following theorem (which is a particular case of a theorem proven 
in [2, p. 495-7]): Let u be harmonic in a bounded domain G, let f be continuous in 
G, let D(u) and D(f) be finite, and let u — f be a B-function. Then u — f ap- 
proaches zero at the boundary if the latter satisfies certain mild conditions. 

Petrovsky [3] has presented a proof that the boundary values are assumed. His 
proof, which also appears in (4, p. 186 ff.], does not make use of (2) at all. How- 
ever, the Petrovsky proof does not appear to extend readily to more general elliptic 
equations, whereas the proof of (2) presented in [1] and the argument presented 
here can be suitably modified so as to apply to other equations. 

Graduate School of Science 


Yeshiva University 
New York 33, N. Y. 


1. R. Courant, K. Frrepricus & H. Levy, “Uber die partiellen Differenzengleichungen 
der mathematischen Physik,’’ Math. Ann., v. 100, 1928, p. 32-74. 
2 


, Courant & D. Hitsert, Methoden der mathematischen Physik, v. 2, J. Springer, 
—_ 987. 


Perrovsky, ‘““New proof of the existence of a solution of Dirichlet’s problem by 
the method of finite differences,’ ’ Uspehit Mat. Nauk., v. 8, 1941, p. 161-170. 

4. J. D. Tamargin & W. FELLer, Partial Differential ’ Equations, Mimeographed lecture 
notes, Brown University, Providence, 1941. 
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1(C, K, L]. Geratp J. Lieperman & Donatp B. Owen, Tables of the Hypergeo- 
metric Probability Distribution, Stanford University Press, California, 1961, 
7 + 726 p., 24 cm. Price $15.00. 


In this volume there are three main tables of the hypergeometric probability 
distribution and a table of logarithms of factorials. The nomenclature of sampling 
inspection is used to describe the parameters of the hypergeometric probability 
distribution. The main tables give the values of p(x) = p(N, n, k, x) and P(x) = 
P(N, n, k, x), where 

N = number of items in a lot, 

n = number of items in a sample taken from the lot, 

k number of defective items in the lot, 

x = number of defective items observed in the sample. 

Then, the probability 


p(x) = Pr{Exactly x defectives are found in sample} 


= ki n! _ (N — k)i (N — 2)! 
(k—z)!(n—az)!z! NUN —k —n +722)!’ 
x being an integer such that (0, n + k — N|) S x S min [n, kj, and P(z) = 
Pr {x defectives or fewer are found in sample} 
= > kin! (N — k)!(N —n)! 
imu (k —a)!(n — alt! NUN —k—n+2)!’ 
where M = max [0,n + k — N}. 

The first table lists the values of p(z) and P(x) to six decimal places for N = 
2(1)49, 50(10)100, n = 1(1) N — 1, k = 1(1) n, ex = O(1) k, for N S 25. 
For N > 25, the values of p(x) and P(z) are given only up to 

N N-1 


ome ¢.2.2* ET: 














N even or odd, respectively. The authors note that by the use of certain symme- 
try relationships, all possible p(z) and P(x) can be obtained. 

The second table gives the values of p(x) and P(x) to six decimal places for 
N = 1000, n = 500, k = 1(1) 500, x = O(1) k/2 (Kk even), (k — 1)/2 (k odd). 
Entries are omitted when p(x) < 10~°. 

The third table gives the values of p(x) and P(x) to six decimal places 
for N = 100(100) 2000,» = 4N,k = n — 1,n, andz = 0(1) n/2 (n even), 
(nm — 1)/2 (n odd). Entries are omitted when p(x) < 10°. 

The fourth table is a list of log N! for N = 1(1) 2000 taken from Logarithms of 
Factorials from 1 to 2000, by D. B. Owen and C. M. Williams, Sandia Corporation 
Monograph SCR-158, December 1959. Values of log N'! in this table are given to 
fifteen decimal places, and were used for the calculation of p(x) and P(x). The values 
of p(x) and P(x) in all tables are claimed to have been computed correct to at 
least eight decimal places before they were rounded to six decimal places. 
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The introductory part of this volume includes the definitions of the hyper- 
geometric function and the various symmetry relationships, applications, approxi- 
mations and interpolations, a summary of some useful formulas on sums of com- 
binatorials, and a bibliography of 66 references. Examples given in applications 
include sequential procedure, test of the equality of two proportions, distribution 
of the number of exceedances, Bayesian prediction, and sampling inspection. 

The reviewer’s immediate reaction to these tables is that the type face is too 
small for easy reading and that the format makes it difficult to find the values of 
the indexing parameters. However, considering the 135,874 entries and the 726 
pages, it would be difficult to eliminate these faults without prohibitive increase in 
both the size and cost of this volume. 


H. H. Kv 


National Bureau of Standards 
Washington, D. C. 


2(G, I, X, Z]. Raupu G. Stanton, Numerical Methods for Science and Engineering, 
Prentice-Hall, Inc., Englewood Cliffs, N. J., 1961, xii + 266 p., 23 em. Price 
$9.00. 


This book is designed as a textbook for an introductory course in numerical 
methods for students in the physical sciences and engineering with a good knowl- 
edge of calculus and differential equations. The selection of topics is fairly standard, 
as one would gather from the following chapter headings: Ordinary Finite Dif- 
ferences, Divided Differences, Central Differences, Inverse Interpolation and the 
Solution of Equations, Computation with Series and Integrals, Numerical Solution 
of Differential Equations, Linear Systems and Matrices, Solution of Linear Equa- 
tions, Difference Equations, Solution of Differential Equations by Difference 
Equation Methods, and the Principles of Automatic Computation. 

The author states that the book was developed from the standpoint of hand 
and desk-calculator techniques, and justifies this on the grounds of his belief that 
“the majority of workers in science and engineering can make great use of numerical 
methods without perhaps ever encountering a problem of sufficient length or 
complexity to justify programming it for an electronic computer.” His final chapter, 
containing only eighteen pages about automatic computation, seems to confirm 
one’s belief that the author views the modern field of numerical computation with 
automatic electronic computers as a spectator rather than as a participant. 


D. M. Y. 


University of Texas 
Austin, Texas 


3[G, S]. Taro Suimpuxu, “General Theory and Numerical Tables of Clebsch- 
Gordan Coefficients,” Progr. Theoret. Phys., Kyoto, Japan, Supplement No. 
13, 1960, p. 1-135. 


General formulas for the Clebsch-Gordan coefficients (j:jsmyme | jijejm), in 
the notation of Condon and Shortley [1], have been given by Wigner and by Racah 


[2], [3]. These formulas are very complex and computationally inconvenient. 
Shimpuku states: “Here we derive a new general expression of C — G coefficients 
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from the theory of spinor representation in three-dimensional rotation group, and 
this expression has a convenient form for practical evaluation (for any given values 
of the parameters).” 

Algebraic formulas for these coefficients, for the special cases j2 = 4, 1, 3, 2 
are given in [1], p. 76-77; similar formulas for 7, = § and 3 are available in sources 
noted in the references in Shimpuku’s paper. Shimpuku tabulates the algebraic 
formulas for j. = 4, 4, $, and 5. 

Numerical tables have been compiled, by Simon at Oak Ridge, for all cases 
where j; S $, je S $. Over 100 pages of numerical tables are given by Shimpuku; 
these tables cover j2 = 5, 4, and 6 for all j, < 6. Each entry is expressed as the 
radical of a rational fraction. 

Shimpuku does not refer to the recent tabulation by Rotenberg et al. [4] of 
3 — j symbols, from which the Clebsch-Gordan coefficients can be readily ob- 
tained. This tabulation covers all values j, S 8, je S 8. However, for the range 
covered by Shimpuku, many users may find his rational fractions more convenient 
than the expressions as products of powers of primes used by Rotenberg. 


GrorGE SHORTLEY 
Booz-Allen Applied Research, Inc. 
Bethesda 14, Maryland 


1. E. U. Connon & G. H. SHortiey, The Theory of Atomic Spectra, Cambridge University 
hae New York, 1935, p. 75. 


E. P. WIGNER, Gruppentheorie und thre Peerage | auf die Quantenmechanik der Atom- 
spekiren, Friedrich Vieweg und Sohn, Braunschweig, 1 Bs 
Racau, “Theory of complex spectra II,’’ Phys. a. v. 62, 1942, p 
: M. RoTENBERG, R. Bivins, N. Merropouis & J. K. Wooren, Jr., The 3-j and 6- 3 Sym- 
bols, Technology Press, Cambridge, 1960. See Math. Comp., v. 14, 1960, p. 382-383, Review 71. 


WA 4{I, X, Z|. Nationa Puysicat Lasoratory, Modern Computing Methods, Second 
Edition, Her Majesty’s Stationery Office, London, 1961, 25 em. Price $3.78. 


This is the second edition of a booklet I praised highly when I reviewed its first 
edition (MT AC, v. 12, 1958, p. 230, Review 96). However, much has changed since 
then, and while I still feel that I shall recommend to every budding numerical 
analyst that he consult this booklet, I must add here to its list of limitations. I 
am tempted to say that it is “modern,” much the same as Gilbert and Sullivan’s 
Major General, but this would be entirely too harsh. 

The booklet contains nothing about linear programming, assignment problems, 
or discrete variable calculations, which play a large role in computation, at least 
in the United States. (Beale in the United Kingdom might claim that these prob- 
lems occur there also.) There is nothing about the Monte Carlo method which is 
very popular, at least in the southwest sections of the United States. (Hammersley 
in the United Kingdom might claim that these problems occur there also.) There 
is essentially nothing (nine lines of text, washing their hands of the whole subject) 
concerning latent roots and characteristic vectors of unsymmetric matrices, al- 
though some of these problems are vital in the study of stability. (This is most 
disappointing of all, for the workers at the National Physical Laboratory were 
spectacular in their early attacks on matrix problems and their reporting of their 
experiences.) There is a tendency to make overly dogmatic statements: “For 
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hyperbolic equations the existence of real and distinct characteristics leads to the 
most satisfactory known method of numerical solution” (p. 105). 

The listing of tables of functions of several variables is not adequately de- 
scribed. In particular, the alluring paper by Kolmogoroff [1] is not mentioned. In 
general, scant attention is paid to important Russian work; the book by Kantoro- 
vich and Krylov [2] is not listed in the bibliography, even though it is available in 
an understandable English translation. 

Despite these criticisms, which might be likened to the disappointment of a 
lover (of the first edition) as his love ages, this is a handy booklet to have available. 
It is a cook book of procedures which are recommended on the experience of a 
perceptive, scholarly, and active computing group. It is less necessary now than 
it was when it was published in its first edition, for SHARE and the other users’ 
groups have made experiences with computers more easily available to other users, 
but this booklet is more precise, less coding-conscious, and more scholarly than 
the reports of the users’ groups. The booklet has been brought up to date on the 
topics it covers; Givens and Householder on latent roots are quoted carefully, in- 
cluding a British interpretation of their impressive work in both avoidance of long 
calculation and analysis of error. On the other hand, many reports of computa- 
tional experience now exist in the literature, which was not the case when the first 
edition was published, so the booklet is no longer a must. 

I would feel unhappy if I knew of this volume and did not have it in my library. 
I suggest that firms which have spent millions of dollars on computers buy a few 
copies even though some isomorph of SHARE is available. 

If a third edition is contemplated, I suggest that the chapter on Finite Difference 
Methods be omitted as non-modern. By implication above, I have suggested 
chapters which should be included. Also a chapter on coding and coding languages 
might reasonably be added. 

I note that there is considerable modernization of outlook (including a chapter 
on Chebyshev series), and this is good. 


C. B. TompxKIns 
University of California 
Los Angeles 24, California 


Se Ko.mocorov, “O predstavlenii nepreryrny kh funktsil neskopkikh peremenny kh 


v vide superpozitsii neprer vaykh funktsil odnogo peremennogo i slozheniia”’ (‘On the repre- 
sentation of a continuous function of several variables in the form of a superposition of con- 
tinuous functions of one variable and their sums’’), Akad. Nauk SSSR, Doklady, 1957, v. 114, 


p. 
2. L. V. Kantorovicn & V. I. KryLov, Approzimate Methods of Higher Analysis, trans- 
lated by Curtis Benster, Noordhoff, Groningen, 1958. 


5(K]. G. W. Rosentuat & J. J. Roppen, Tables of the Integral of the Elliptical 
Bivariate Normal Distribution over Offset Circles, LMSD-800619, Lockheed 
Missiles and Space Division, Sunnyvale, California, May 1961, iii + 92 p., 28 
em. 


These tables give the probabilities of being inside various circles not about the 
mean from a bivariate normal distribution having unequal variances. The range 
of the tables includes values of the mean up to three times the standard deviation. 


AUTHOR’s SUMMARY 


a mem 
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6{K, X]. V. V. SoLopovnixov, Introduction to the Statistical Dynamics of Automatic 
Control Systems, Translation edited by John B. Thomas and Lotfi A. Zadeh, 
Dover Publications, New York, 1960, xx + 307 p., 20 cm. Price $2.25 (Paper- 
bound). 


This book, first published in Russian in 1952, gives an excellently written, self- 
contained account of the principles of the analysis of linear systems, the statistics 
of random signals, and the theory of linear prediction and filtering. The translation 
is well done. In addition to treating exact methods, the author discusses methods 
of obtaining approximate solutions to various problems. 

The first three chapters are devoted to a discussion of the theory of the transients 
in a linear system produced by deterministic signals, to the elements of probability 
theory, and to the basic concepts of the theory of stationary random processes. 

Chapter IV discusses the criterion of least mean-square error. Linear and square- 
law detectors are used to show how some nonlinear systems may be treated. 

In Chapter V the problem of using numerical methods to approximate spectral 
distribution curves is treated. 

Chapters VI, VII, and VIII contain the derivation and application of formulas 
from which one may obtain the transfer function yielding a minimum mean-square 
error from the knowledge of the spectral densities of the signal and noise. The last 
of these chapters treats the case where the signal is composed of two parts, one 
deterministic and one random. 

The book contains four appendices. Appendix I consists of five-place tables of 


the functions = and = * for x = 0(.01)10.0(.1)20(1)100. Appendix II contains 


tables of the first five Laguerre functions to five significant figures for values of the 
argument in the range 0(.01)1.0(.1)20(1)30. Appendices IIIa and IIIb give five- 
place tables for the calculation of the so-called phase characteristic function from 
straight-line approximations of the logarithm of the spectral-density function. 
Appendix IV gives a table of integrals 


bes ate olde getone 
be = 355 | ape)" ~ ONT 








where 
Gn(jw) = bo(jw)” + bi(jw)”* + --- + dn, 
H,(jw) - Ao(jw)" + Ai(jw)”"" + sate + A, ? 
and all roots of H,(jw) are in the upper half-plane. 
A. H. T. 


VY 7[L]. Davin J. BenDantet & Wiriu1am E. Carr, Tables of Solutions of Legendre’s 


Equation for Indices of Nonintegral Order, University of California Lawrence 
Radiation Laboratory, Livermore, UCRL-5859, September, 1960, 68 p., 28 cm. 
Available from the Office of Technical Services, Washington 25, D. C. Price 
$1.75. 


We employ the usual notation for hypergeometric and Legendre functions [1]. 
Let 


f(x) = oFy(—v/2, v/2 + 4;4;2°); fo(x) = mF i(4 — v/2, 1 + v/2; 3/2; 2°). 








118 REVIEWS AND DESCRIPTIONS OF TABLES AND BOOKS 


Then f(z) and f2() satisfy the differential equation 
(2° — 1) d’f/ de® + 2x df/ dx — v(v + 1)f = 0, 
and we have 
Px) = arfi(z) + aefo(x), Q(x) = bifi(x) + beafe(z), 
where 
a, = we (P($ — »/2)T(1 + v/2)]", as = —2e'?[P(4 + »/2)P(—»/2)]7 


ai ietss he ’T(a + v/2) sin vx/2 “— eri + v/2) cos vr/2 
7 T(1 + v/2) Bie T'(3 + v/2) 
Table 1 gives fi(z) to 5S, corresponding to x = 0(0.01)0.99, » = 0(0.0625) 
1(0.125)10(0.25)36. 


For a given y, let a” be the i-th zero of f;(x). Then Table 2 gives 4S values of 
a” and of the integrals 








[ fi(x) dz, £" {fi(x)}* dz, for v = 0.4375(0.0625)36. 
0 0 


Tables 3 and 4 present corresponding data for f2(z). 

The tables are essentially new, and were obtained on an automatic computer. 
An introduction gives a few definitions (the coefficients a, and a2 in the formula 
for P,(x) contain typographical errors). There is no discussion of formulas used 
to perform and check the calculations. No attempt is made to give closed-form 


results, which would be useful to the applied worker. For example, we can show 
that 


fiz) = (cos 6/2)~™* cos ims + say (tan 0/2 + 0/2)\ {1 + O(1/s*)}, 
s¥1 
| fi(t) dt = (2N,)~*(cos 0/2)" sin {No a vay. (—3 tan 6/2 + 0/2) 
0 LV 
-{1 + O(1/r’)}, 


fo(x) = (2N2)~*(cos 6/2) *” sin {Neo + TaN: (tan 6/2 + 09/2) 


-{1+ O(1/**)}, 
/ f2(t) dt = —(2N2)~*(cos 0/2)"” cos {Neo + a (—tan 6/2 + 36 a} 
0 LVe 
-{1 + O(1/*)} 


NY = o(v + 1)/4, No = (v —1)(v + 2)/4, x = sin0/2, O< O0<-. 


Indeed, using these results, we have spot checked the entries. Tables 1 and 3 and 
values of the zeros of f; and f. in Tables 2 and 4 appear to be correct. 


If vy = 20, and a = 1, 2, 3, the values of | fi(t) dt are erroneous, as is also 
0 


Z1 
the value of [ fi'(t) dt. It is curious that the error in each case is about 0.0021. 
0 





a wm aaa ae ee 


ea = 





REVIEWS AND DESCRIPTIONS OF TABLES AND BOOKS 119 


Similar check calculations for f.(z) reveal a persistent error of about 0.6024. Thus 
the tables of the numerical values of the integrals should be used, if at all, with 
caution. We have corresponded with one of the authors (D.J.B.). He has checked 
those entries in Tables 1 and 3 against known values of Legendre polynomials and 
finds that they are correct. The reason for the bias in the values of the integrals is 
not known, but he suspects that it arises from the binary-to-decimal conversion. 
We conclude with the “trite” observation that aut~—atic computers cannot be 
trusted implicitly, and that the need for analysis and checking remains. 


Y. L. L. 


Midwest Research Institute 
Kansas City, Missouri 


1. A. Erp&.y1, et al., Higher Transcendental Functions, Vol. 1, McGraw-Hill, New York, 
1953 . 


8[L}. Lupo K. Freve, & J. W. Turtey, Tables of Iterated Sine-Iniegral, The 
Dow Chemical Company, Midland, Michigan, 1961. Deposited in UMT File. 


Three tables of decimal values of the iterated sine-integral, Si (x), are herein 
presented, as computed on a Burroughs 220 system supplemented by Cardatron 


equipment, which permitted on-line printing of the results in the desired tabular 
format. 


Table 1 presents the values of Si (x) to 9D for n = 1(1)10, x = 0(0.2)10. 
Table 2 gives values of this function to 7D for n = 0(0.05)10, x, 2x, 3x, and 
Table 3 gives for n = 1(1)10 the values to 9D of the first thirty extrema, which 
correspond to x = mz, where m = 1(1)30. 

In an accompanying text of three pages the authors describe in detail the 
method of calculation and the underlying mathematical formulas. It is there 
stated that the entries in Table 2 were computed to 9D prior to rounding. The 
entries in Table 3 are claimed to be accurate to within a unit in the final decimal 
place, and the authors imply in their explanatory text that comparable accuracy 
was attained in the computation of the entries in Table 1. 

The tabular data corresponding to the values of n different from unity constitute 
an original contribution to the literature of mathematical tables. 


J.W. W. 


O[L, X]. Hans Sacan, Boundary and Eigenvalue Problems in Mathematical Physics, 
John Wiley & Sons, Inc., New York, 1961, xviii + 381 p., 24 em. Price $9.50. 


This attractive newcomer to the ranks of the textbooks on methods of mathe- 
matical physics comes to us directly from Moscow (where, for the past four years, 
the author has been an Associate Professor of Mathematics at the University of 
Idaho). This book contains material which has been used in the author’s classes 
to seniors and beginning graduate students in mathematics, applied mathematics, 
physics, and engineering for the past five years. The author’s stated purpose is not 
to present a vast number of seemingly unrelated mathematical techniques and 
tricks that are used in the mathematical treatment of problems which arise in 
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physics and engineering, but rather to develop the material from a few basic con- 
cepts; namely, Hamilton’s principle together with the theory of the first variation, 
and Bernoulli’s separation method for the solution of linear homogeneous partial 
differential equations. The author’s persuasive style appears certain to gain ad- 
herents for his viewpoints on many college campuses this coming fall. 

Hamilton’s principle and the theory of the first variation occupy Chapter 1. 
The representation of some physical phenomena by partial differential equations 
(vibrating string and membrane, heat conduction and potential equation) forms 
the subject matter of Chapter II. Chapter III contains general remarks on the 
existence and uniqueness of solutions and the presentation of Bernoulli’s method 
of separation of variables, while Chapter IV is devoted to Fourier series. Chapter 
V deals with self-adjoint boundary-value problems, the concept of their eigenvalues 
being developed according to the elementary method of H. Pruefer in Mathe- 
matische Annalen, v. 95 (1926). Chapters VI and VIII, on special functions, deal 
with Legendre polynomials and Bessel functions, and spherical harmonics, respec- 
tively. Chapter VII develops the characterization of eigenvalues by a variational 
principle; while the final Chapter [X is devoted to the nonhomogeneous boundary- 
value problem (Green’s function and generalized Green’s function). 

The text is well designed for class room use. The author intends it to be used 
in a two-semester three-credit course. Each chapter is generously provided with 
interesting exercises (answers and hints are provided at the end of the book for the 
even-numbered problems). A recommended supplementary reading list concludes 
each chapter. A welcome innovation is the detailed appendix, containing a con- 
densation of topics with which “the student who wishes to take this course with a 
reasonable chance to succeed should be familiar.” 


J. B. Draz 
The University of Maryland 


College Park, Maryland 


10(S]. Caartes DeWitt Coteman, WituiaM R. Bozman & WiiuiaM F. MEaGERs, 
Table of Wavenumbers, Volumes I and II, U. S. Department of Commerce. 
Volume I—2000 A to 7000 A, and Volume II—7000 A to 1000 u, 1960, vii + 
500 p., and vii + 534 p., 35 em. Price $6.00. 


A two-volume table for converting wave lengths in standard air to wave num- 
bers in vacuum was computed by using the equation ovae = 1/(mdAgir), where n 
was computed from Edlen’s 1953 equation for the refractive index of air. Wave 
numbers are given to the nearest 0.001 K (cm™") for wave lengths from 2000 to 
7000 A in volume I, and 7000 A to 1000 u in volume II. Proportional tables are 
given for linear interpolation between entries of \. Also included are the vacuum 
increase in wave length, (n — 1), and the refractivity of standard air in the form 
(n — 1) X 1000. 


AvutTHors’ SUMMARY 


11[W]. Guy H. Orcutt, Martin GREENBERGER, JOHN Korpet & Auice M. 
Rivurn, Microanalysis of Socioeconomic Systems: A Simulation Study, Harper & 
Brothers, New York, 1961, xviii + 425 p., 21 cm. Price $8.00. 


In this book the authors discuss an experimental calculation carried out on a 
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high-speed calculator for the purpose of predicting the population trend during 
the period 1950 to 1960. The calculation is based on an elaborate model which is 
designed to simulate the demographic characteristics of the population by means 
of a large number of typical household units (approximately 5000). Each household 
unit represents a segment of the population, such as the married white female 
members between the ages of 20 and 25. The calculation proceeds in short time incre- 
ments (months) until the final state is reached. The distribution of the population 
is computed at each time interval, taking into consideration the probabilities for 
such occurrences as births, deaths, divorces, etc. This process may be compared to 
the use of the Monte Carlo method for the solution of the diffusion equations in 
physics. 

The reviewer believes that the authors would have better served the interests 
of future research in this field if they had devoted their discussion to a factual 
description of the results attained and difficulties encountered in carrying out this 
interesting but rather restricted experiment. However, as indicated by the some- 
what pretentious title, this is not the primary emphasis of the book. The authors 
appear to stress the potential application of their techniques in the simulation of 
the total social-economic structure of the United States; and the book is promoted 
as a “pioneer work with a new approach to the scientific study and analysis of 
social systems, employing the major tools of modern research.”’ The enthusiasm 
of the authors for their methods would have been more easily understandable if 
their calculations would have accurately predicted what the population distribu- 
tion will be in 1970, rather than what it was in 1960. 

| A 


12(W, X]. Metvin Dresuer, Games of Strategy: Theory and Applications, Prentice- 
Hall, Inc., Englewood Cliffs, New Jersey. 1961, xii + 186 p., 23 cm. Price 
$9.00. 


This small volume on zero-sum two-person games contains essentially the whole 
story on finite games and a great deal on infinite games. It can be profitably read 
by anyone with some calculus and the first chapter or so of matrix theory behind 
him. The author presents an elementary proof of the minimax theorem which also 
yields a good computational procedure for solving finite games. The properties of 
optimal strategies are then discussed in an exhaustive and illuminating manner, 
and various methods of solving games are described. The subject of infinite games, 
filling one-half the book, is treated next, and the topics covered include games 
with convex payoff functions, games of timing, and games with separable payoff 
functions. Numerous examples of such games, described in military terms, are 
given and their solutions discussed thoroughly. The author’s style is pleasant, and 
the printing and layout of the book are attractive. 


JosEPH BRAM 
Applied Mathematics Laboratory 
David Taylor Model Basin 
Washington 7, D. C. 


13[X, Z]. L. Layton, H. Smita & L. Cuatrieip, Proceedings of Executive Seminar 
on the Use of High-Speed Calculators for the Solution of Naval Problems, Applied 
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Mathematics Laboratory, David Taylor Model Basin, Washington 7, D. C., 

DTMB Report 1519, May 1961, iv + 355 p., 27 cm. 

This book consists of twenty unclassified papers presented at a seminar at the 
Applied Mathematics Laboratory, David Taylor Model Basin, Carderock, Mary- 
land, during 7-9 September 1960. Six additional papers classified ‘“‘Secret’”’ and one 
classified “Confidential” are not included in this volume. 

The papers are oriented toward the use of high-speed computers in the solution 
of Naval Problems, with emphasis on applications drawn primarily from the Bureau 
of Ships activities. The general areas covered are: (1) engineering research, (2) 
management data analysis, (3) large-scale data processing, (4) operations research, 
and (5) tactical and strategic planning. 

The text is double-spaced and easy to read; however, the quality of reproduc- 
tion of the photographs leaves much to be desired. 

There is not enough space to review each paper separately, so that the following 
statements may do some injustice to individual papers. Several authors report on 
their own practical experience, and do not give a perspective to the subject dis- 
cussed. However, there are many excellent papers, especially “(Computer Tech- 
nology Outside the USA” by Dr. 8. N. Alexander, “Nuclear Reactor Design Cal- 
culations” by Joanna Wood Schot, ““Mathematical Calculation of Shiplines” by 
Dr. F. Theilheimer, ‘“The Solution of Naval Problems on High-Speed Calculators’’ 
by Dr. H. Polachek, and several others. The paper, ‘“On Teaching of Mathematics,” 
by Dr. Francis D. Murnaghan, should be read by every mathematics teacher. All 
in all, the book offers valuable reading for both the beginner and the experienced 
computer specialist. 

It is unfortunate that the remarks of the keynote speaker, Professor Howard 
H. Aiken were not recorded for this volume, since he is recognized as the father of 
modern computers. 


RatpH A. NIEMANN 
U. S. Naval Weapons Laboratory 


Dahlgren, Virginia 


14[Z]. WitHELM KAMMERER, Ziffernrechenautomaten, Akademie-Verlag, Berlin, 

1960, viii + 303 p., 24 em. Price DM 29. 

This well-written book is based on a course of lectures given by the author at 
the Friedrich-Schiller University, in Jena. It discusses computer components, their 
organization into the various organs of a computer, the logical organization of 
computers, and the fundamentals of programming. 

The first chapter discusses the binary number system and Boolean algebra. 
The second chapter is concerned with the nature of arithmetic operations and 
methods for realizing them by automatic devices. 

Chapter three deals with the structure of an automatic computer and the re- 
quirements that must be imposed on it. It illustrates these requirements and 
methods by which they have been satisfied, by reference to various computers. 

Chapter four discusses in some detail various well-known computer components 
and methods for organizing them into computer organs. The newer components 
are not treated. 
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The final chapter, Chapter five, is devoted to the principles of programming. 
Several illustrative problems are coded for an imaginary single-address machine. 
The problems of using a library of subroutines are discussed, as are the notions of 
relative addresses and floating addresses 

The bibliography given at the end of the book is not an extensive one. The oldest 
references in it are dated 1951. This is somewhat unfortunate, for the reader cannot 
gain any impression therefrom of the historical development of the subject. The 
omission of any reference to the fundamental work of von Neumann on computers 
is, to the reviewer, a great oversight. 

A. H. T. 


15|Z]. Herpert D. Leeps & Geratp M. WEINBERG, Computer Programming 
Fundamentals, McGraw-Hill Book Company Inc., New York, 1961, ix + 368 
p., 23 em. Price $8.50. 


Nominally an introductory textbook on digital computing techniques and 
applications, this book presents a readable account of the basic principles of pro- 
gramming and coding for a specific machine, namely, the IBM 7090 computer. 
No mathematical knowledge beyond elementary algebra is required. The first 
section delineates the fundamental characteristics and special capabilities of a 
computer and then highlights the preparatory steps required to obtain a machine 
solution. The longer second section is devoted to an exposition of flow-diagram- 
ming and coding for the IBM 7090 computer. 

In view of the fact that the book is addressed to “students in business adminis- 
tration, economics, and other nontechnical fields as well as the physical sciences 
and mathematics courses’, the authors are disappointingly vague on the subject 
of programming techniques and procedures for the solution of large-scale data 
processing problems. Such significant developments as business compilers (COBOL, 
IBM Commercial Translator, etc.), sort generators, and report generators are not 
even mentioned. The value of the book as a general text on computer fundamentals 
is further lessened by the omission of references and supplemental readings. Con- 
sequently, the reviewer believes that this volume will be primarily suitable as a 
general IBM 7090 programming manual for nontechnical readers. It is written in 
a lively, lucid style that can be easily comprehended by the layman. 


Mitton SIEGEL 
Applied Mathematics Laboratory 
David Taylor Model Basin 
Washington 7, D.S. 


16[Z]. W. W. Pererson, Error-Correcting Codes, The Technology Press and John 
Wiley & Sons, Inc., New York, 1961, x + 285 p., 24 em. Price $7.75. 


The journal literature on algebraic coding theory has become so extensive 
lately that a book has been needed to give perspective and order to the field. This 
excellent book not only fills this need but also improves greatly on the presentation 
in many of the journal articles. In conjunction with the literature on probabilistic 
schemes of coding and decoding, Peterson’s book gives an essentially complete 
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picture of coding theory as currently known. Most of the material does not appear 
elsewhere in book form, and a considerable amount is original. 

The style of writing is remarkably successful in developing insight and intuition 
without appreciably sacrificing rigor or conciseness. The book is almost self-con- 
tained and includes a development of the required algebraic concepts and theorems. 
A discussion of ways to implement algebraic operations, particularly on poly- 
nomials and Galois field elements by shift register circuits, should help the engineer 
to understand and use modern algebra, both in coding theory and elsewhere. 

The first part of the book discusses linear codes, which are group codes or 
parity-check codes generalized to include non-binary alphabets. This includes a 
general treatment, some theoretical bounds on error-correcting ability, and a dis- 
cussion of several specific classes of linear codes. Next, after some mathematical 
development, the theory and implementation of cyclic codes is discussed from 
several different viewpoints. 

Bose-Chaudhuri codes, which are the most important of the known algebraic 
codes, are elegantly treated after cyclic codes. A simple derivation of their error- 
correcting capabilities is given, and two decoding techniques are presented. The 
remainder of the book treats burst-error correction, other approaches to decoding, 
recurrent codes, and the checking of arithmetic operations. 

The appendices include a table of irreducible polynomials over the field of two 
elements. They are arranged in order as minimum polynomials of the elements of 
Galois field, and this makes it possible to find generator polynomials for Bose- 
Chaudhuri codes almost by inspection. 

The book is highly recommended to engineers and mathematicians interested 
in coding, information theory, communication, and computers. 


R. G. GALLAGER 
Massachusetts Institute of Technology 


Cambridge 39, Massachusetts 


17(Z]. A. Unear, Proceedings Editor, Proceedings of the 1959 Computer Applica- 
tions Symposium, sponsored by the Armour Research Foundation of Illinois 
Institute of Technology, Chicago, 1960, x + 155 p., 23 em. Price $3.00. 


Like most symposia, this little book contains quite a variety of papers, some of 
them excellent, and some of them only mediocre. The paper entitled “Fortran 
Experience and Remote Operation by Non-computer Specialists” in conjunction 
with its subsequent panel discussion is, in the opinion of this reviewer, alone well 
worth the price of the book. This is by no means the only interesting paper. Each 
paper is followed by the type of bantering discussion that usually takes place in a 
meeting of a group of specialists. 

The papers include: “Shareholder Record-Handling with the Aid of Character- 
Recognition Equipment,” “Around the World in Eighty Columns,” “Cost Reduc- 
tion Through Integrated Data-Processing,” “Some Aspects of Computer Tech- 
nology in the USSR,” “Experience and Plans for Marketing-Research Operations,” 
“A Modern Approach to Inventory Control Utilizing a Large-Scale EDPM,” 
“Current Developments in Common-Language Programming for Business. Data 
Systems,” “Linear Programming on the Bendix G-15 Computer,” “The Design 
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and Use of the APT Language for Automatic Programming of Numerically Con- 
trolled Machine Tools,” “A Quasi-Simplex Method for Designing Suboptimum 
Packages of Electronic Building Blocks,” “The International Algebraic Language 
and the Future of Programming,” “Training for Engineers and Scientific Applica- 
tions via Compilers, Interpreters, and Assemblers,” “Scientific Design Procedure 
Utilizing a Small Computer,” and “FORTRAN Experience and Remote Operation 
by Non-Computer Specialists.” 

This inexpensive book is recommended for inclusion in the libraries of compu- 
tation laboratories and of individual programmers. 


Ricuarp V. ANDREE 
University of Oklahoma 
Norman, Oklahoma 





TABLE ERRATA 


307.—E. P. Apams & R. L. Hippistey, Smithsonian Mathematical Formulae and 
Tables of Elliptic Functions, second and third reprints, The Smithsonian In- 
stitution, Washington D. C., 1947 and 1957. See also MT AC, v. 12, 1958, p. 262 
and earlier references cited there. 


The following corrections should be made in section 6 ‘“0, on p. 126: 


9.6 9.6 
Formula 4: In the series for e”” *, for + = , read — > 
: >: 


4 4 
Formula 8: In the series fore" *, for + or wadlo oe 


24° 

The corrected series agree with those appearing in B. O. Peirce, A Short Table 
of Integrals, third revised edition, Ginn & Co. Boston, 1929, p. 92-93. 
T. H. SourHarp 


Alameda State College 
Hayward, California 


Ep1tToriat Norte: For further information regarding the expansion of en | +, see Richard 
Kelisky, ‘“The numbers generated by exp (arctan z),” Duke Math. Journal, v. 26, 1959, p. 
569-581, wherein the coefficients of the first twelve powers of z in the expansion are given. 


CORRIGENDA 


M. Ascuer, “Explicit solutions of the one-dimensional heat equation for a 
composite wall,”’ Math. Comp., v. 14, 1960, p. 346-353. 


The last sentence of the derivation of the estimate of error (on page 350) yields 
the result O( Ar) rather than O( Az’). Nevertheless, the error is O( Az’), as shown 
by the numerical example and by abstract No. 576-136, Amer. Math. Soc. Notices, 
v. 7, no. 7, December 1960, p. 944. 


MartTIN GREENBERGER, “An a priori determination of serial correlation in com- 
puter generated random numbers,” Math. Comp., v. 15, 1961, p. 383-389. 


On p. 387, the first two sentences in the paragraph immediately following equa- 
tion (19) should follow instead equation (11) on p. 385. The reference to equation 
(1) in the first of these sentences should be changed instead to equation (10). 


Sipney Kravitz, “Divisors of Mersenne Numbers 10,000 < p < 15,000,” 
Math. Comp., v. 15, 1961, p. 292. 


In the 7th column of the table, 8th line down, the number 12,063 should read 
13,063. 
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